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SECTION  I 


INTRODUCTION 

The  development  of  transonic  prediction  methods  to  complex  wing-fuselage-pylon- 
stores  combinations  is  of  great  importance  because  of  its  potential  applications  in  fighter 
aircraft  design,  and  for  the  detailed  flow-field  analysis  required  to  ascertain  the  safe 
release  of  and  the  subsequent  accurate  prediction  of  the  trajectory  of  the  external  stores 
after  release.  Several  finite-difference  and  finite-volume  computer  codes  developed  by 
Jameson  and  Caughev  (Reference  1-3)  viz.,  FLO-27,28, and  30  can  predict  transonic 
flow  including  the  approximation  of  embedded  shock  wave  location  and  strength,  about 
a  swept  wing  attached  to  either  a  wall  (plane  of  symmetry)  or  a  simplified  fuselage 
geometry.  The  simplified  fuelage  geometry  varies  from  an  infinite  cylinder  to  a  finite 
cylinder  of  moderately  varying  cross  sections  in  these  codes.  The  extension  of  the 
procedure  for  wing-fuselage  treatment  in  these  codes  to  the  wing-fuselage-pylon-store 
case  is  not  straightforward  due  to  several  reasons.  Recent  numerical  experiments  in  this 
direction  to  investigate  the  direct  approaches  have  not  been  very  successful. 

Strong  viscous-inviscid  interactions  in  transonic  flow  regime  introduce  additional 
complications  to  the  exact  prediction  of  the  fluid  motion.  The  existence  of  shock  and  the 
shock-induced  separation  regions  in  transonic  flows  is  quite  common.  Although  precise 
mathematical  procedures  for  the  solution  of  the  complete  set  of  equations  describing  the 
fluid  flow  are  available  in  the  literature,  they  are  all  concentrated  to  isolated  individual 
areas.  In  the  absence  of  a  comprehensive  numerical  technique  for  a  complex  problem, 
attempts  to  include  appropriate  elaborate  solutions  locally,  along  with  the  global  ap¬ 
proximate  solutions,  have  been  very  rewarding.  Several  of  these  ideas  have  been  in 
practice  in  the  industry,  A  very  good  example  for  the  above  is  the  McDonnell  Douglas 
Aircraft  Company's  (MCAIR)  Equivalent  Simple  Body  (ESB)  concept  as  discussed  in 
Reference  4. 

The  validity  of  these  approximations  depend  upon  the  significance  of  each  of  the 
individual  areas.  A  first  estimate  of  their  contributions  to  the  solution  could  be  obtained 
by  the  asymptotic  or  matched  asymptotic  expansion  principles  with  single  or  multiple 


parameters  representing  the  complexity  of  the  problem.  Earlier  efforts  in  this  area  to 
simulate  a  complex  flow  field  with  Transonic  Small  Disturbance  equations  proved  to 
be  very  limited  (Reference  5).  With  the  above  classifications  in  mind,  we  can  observe 
that  the  exact  and  extensive  numerical  computation  of  a  certain  simplified  equation  of 
motion  for  the  transonic  flow,  that  has  significant  contribution  from  at  least  more  than 
two  categories,  may  not  yield  accurate  answers. 

In  the  present  report,  a  methodology  to  solve  the  above  problems  is  intended.  An 
integral  equation  method  for  the  solution  of  the  transonic  small  disturbance  equation  is 
applied.  The  method  involves  a  non-linear  correction  term  applied  to  the  basic  linear 
solution  obtained  with  the  suitable  compressible  coordinate  transformation.  The  solution 
process  as  discussed  in  Reference  6  is  also  highlighted  in  Section  II.  Viscous  effects  are 
included  by  calculating  the  two-dimensional  compressible  boundary  layer  displacement 
thickness.  A  displaced  body  concept  or  Lighthill's  equivalent  blowing  concept  could 
be  used  to  incorporate  the  viscous  effects  both  on  the  non-lifting  and  lifting  bodies. 
The  wake  effects  could  also  be  introduced  in  a  similar  way.  Methods  of  introducing 
the  separation  effects  are  discussed,  but  the  modification  of  the  code  to  incorporate  the 
separation  effect  is  not  made.  A  brief  review  of  the  basic  scheme  for  viscous-inviscid 
interaction  is  outlined  in  Section  III. 

An  approximate  method  for  the  flow-field  velocity  computation  for  the  separated 
store  arrangement  has  been  developed.  The  experimental  data  for  configurations  101 
and  102  (see  Appendix  A)  showed  regions  of  shock  and  shock-induced  separation.  In 
the  present  method  the  shock  location  has  been  computed  by  a  forward  and  backward 
matching  procedure  (in  the  subsonic  regimes)  in  the  solution  process  of  two-dimensional 
Navier-Stokes  equations.  A  predictor-corrector  MacCormack’s  explicit  scheme  has  been 
used  to  compute  the  two-dimensional  internal  flow  between  the  wing  and  the  store. 
The  location  and  strength  of  the  shock  so  computed  show  good  agreement  with  the 
experimental  values.  The  solution  methodology  along  with  the  numerical  procedure  is 
briefed  in  Section  IV. 

Computation  of  the  store  pressure  is  another  aspect  discussed  in  this  report.  Again, 
the  computation  of  on-body  pressures  for  complex  configurations  is  a  challenging  task. 


I 


The  rapid  decay  of  disturbances  that  could  be  expected  at  the  off-body  points  and  hence 
the  validity  of  making  a  lower  order  approximations  of  the  full  Gas  Dynamic  equations 
are  no  longer  valid  for  the  on-body  computations.  But,  due  to  the  low  cross  flow  ' 
components  of  the  flow  variables  at  low  angles  of  attack,  the  axi-symmetric  nature  of 
the  flow  (assuming  that  the  store  is  axi-symmetric)  around  the  store  relaxes  the  transonic 
non-linear  effects  appreciably.  Hence,  based  on  these  experiences  and  due  to  the  fact  that 
there  is  no  comprehensive  and  accurate  computational  procedure  available  for  complex 
geometries  to  predict  the  on-body  pressures,  yet  another  semi-rigorous  procedure  has 
been  developed  for  the  store  pressure  calculations.  The  proposed  procedure  has  been 
discussed  in  some  detail  in  Section  V. 


SECTION  n 


INTEGRAL  EQUATION  METHOD 

The  three-dimensional  Transonic  Small  Disturbance  equation  for  the  velocity  poten¬ 
tial,  <fi,  can  be  written  as: 

(1  ~  +  <t>yy  4-  <t>zz  =  K<t>x<t>xx  (1) 


where 
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and  the  subscripts  x,  y,  and  z  denote  the  differentiation  with  respect  to  the  respective 
variables,  7  is  the  ratio  of  specific  heats,  and  U0 0  ,  are  the  free  stream  velocity  and 
Mach  number,  respectively. 

Using  Green’s  Theorem  and  applying  Goethert’s  Transformation  and  after 
simplifications  (Reference  6),  the  velocity  components  in  the  three  different  directions 
are  written  as: 
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and  G^,  Vl,  U'l  are  the  linear  velocities  in  the  S,  and  2  directions. 

Again,  after  introducing  a  transformation  and  proceeding  as  outlined  in  (Reference  5) 

$  =  Tan  (|ti)  ;  r)  =  Tan  (|t2)  :  ?  =  Tan  (Jf8)  (5) 

we  obtain,  a  quadratic  equation  in  G, 

Z  =  tiL+  y  -  y  (6) 

where  7a  represents  the  triple  integral  term  in  equation  (2). 

A  6-point  Gauss  quadrature  method  is  used  to  numerically  integrate  /n.  This  gives 
the  solution  of  the  non-linear  0  velocity, 

a  =  1  -  - W  - /«)  (7) 

for  local  subsonic  flow,  and 

u  =  1  +  (8) 


for  supersonic  local  flow. 

The  sonic  point  is  determined  by  the  vanishing  of  the  descriminant  and  the  smooth 
transition  of  from  subsonic  to  supersoinc  flow  criterion.  No  predetermined  position  of 
the  shock  is  necessary  for  the  non-linear  flow  computation.  The  calculation  of  G  and  H 
velocities  involves  only  the  evaluation  of  the  /„  and  /*  integrals. 

G(2,  V,  2)  =  Gl  -  h  (9) 


and 

Q}($,tl,2)  =  wl  —  /® 


(10) 


The  computation  of  the  linear  velocities  is  done  by  the  Douglas  Neuman  type  panel 
method  (Reference  7).  The  entire  body  is  represented  by  linear  surface  source  and  doublet 


singularities  to  effectively  simulate  the  thickness  and  lifting  effects, respectively  (Figure 
1).  A  higher  order  panel  incorporating  a  parabolic  vorticity  distribution  could  also  be 
incorporated  (the  present  code  has  this  capability).  The  inclusion  of  these  higher  order 
panels  becomes  effective  only  for  thick  wings  with  fairly  large  leading  edge  radii. 

The  panel  method,  although  cumbersome  since  it  requires  the  entire  body  to  be 
represented  by  quadrilateral  panels,  has  several  desirable  features.  The  present  store 
pressure  computation  procedure  follows  the  steps  mentioned  below: 

(a)  It  is  an  integral  method  of  solution  of  the  Laplace  Equation  and  hence  fully  con¬ 
servative. 

(b)  Suitable  for  arbitrary  geometry  of  the  body. 

(c)  Lifting  and  non-lifting  bodies  could  be  handled  simultaneously  with  as  many  num¬ 
bers  of  lifting  and  non-lifting  sections  as  desired. 

An  important  point  to  be  observed  in  panelling  technique  is  that  whenever  there  are 
multiply  connected  regions  in  the  domain  of  computation,  it  is  necessary  to  provide  an 
artificial  cut  at  an  appropriate  location  to  make  the  region  simply -connected.  Artificial 
cuts  should  be  made  at  appropriate  locations  that  shall  not  disturb  the  flow  appreciably. 
This  is  necessitated  due  to  the  fact  that  the  application  of  Green’s  theorem  during 
the  derivation  of  the  integral  equation  to  obtain  the  strength  of  source  and  doublet 
singularities  requires  the  region  to  be  simply-connected.  This  is  a  definite  experience 
gained  during  the  course  of  the  present  investigation.  Examination  of  the  integrals 
In.  h,  I  a  would  show  that  the  non-linear  contribution  is  significant  only  near  regions  of 
large  gradients  of  the  linear  velocities,  <pxx.  Hence,  at  the  regions  of  rapid  accelerations, 
the  non-linear  method  predicts  higher  shock  strengths  compared  to  the  physically  correct 
values.  This  was  observed  in  the  cases  of  configurations  101  and  102.  A  Fortran  program 
listing  of  the  non-linear  integral  equation  correction  is  given  in  Appendix  B. 

As  the  transonic  parameter 
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become  larger,  particularly  when  the  Mach  number  is  increased  with  the  thickness  r 
being  fixed  (e.g.,  configurations  101  and  121  at  Moo  =  0.975  ),  the  linear  velocities 
considerably  overpredict  the  experimental  values.  This  is  due  to  excessive  stretching  of 
the  co-ordinates  when  the  compressibility  corrections  are  made  to  the  linear  velocity. 
One  way  to  overcome  this  difficulty  is  to  reduce  the  free  stream  Mach  number  in  the 
linear  calculation.  An  alternative  method  is  to  introduce  an  artificial  term  in  the  non¬ 
linear  correction  procedure  that  will  reduce  this  effect.  An  attempt  in  this  direction  was 
made  by  the  principle  of  asymptotic  expansion  with  a  small  parameter  involving  &.  No 
correct  estimate  could  be  derived  for  various  Mach  numbers  and  thicknesses.  Hence,  all 
the  linear  velocity  calculations  were  done  at  the  same  transonic  Mach  numbers  as  the 
experimental  values.  Section  IV  describes  a  local  internal  flow  calculation  procedure  for 
fixing  this  problem. 

t 

It  is  interesting  to  note  that  in  those  configurations  where  the  wing  has  two  degrees 
angle  of  attack  (configuration  201,  221,  etc.)  the  computed  results  show  overall  good 
agreement  with  the  experiment.  But  with  the  0  degree  angle  of  attack  of  the  wing,  regions 
of  shock  are  noted.  This  is  also  shown  in  the  shadowgraph  experiments  of  Reference  8. 
The  appearance  of  shock  is  due  to  the  increased  acceleretion  on  the  bottom  surface  of 
the  wing  (e.g.,  101,  102,  etc.).  At  Moo-  =  1-025,  the  linear  solutions  were  obtained  by 
using  a  free  stream  Mach  number  of  0.9522,  which  corresponds  to  the  downstream  Mach 
Number  across  a  normal  shock.  It  is  indeed  questionable,  but  the  results  thus  obtained 
show  fairly  good  agreement  since  the  detached  shock  is  sufficiently  far  upstream  and 
is  nearly  normal  and  hence  the  entire  panels  are  within  the  subsonic  zone.  At  higher 
supersonic  Mach  numbers,  however,  these  problems  would  be  much  more  complicated 
not  only  because  of  the  interaction  of  shocks  emanating  from  the  wing  leading  edge  and 
the  store  nose,  but  also  of  the  presence  of  panels  that  are  lying  outside  the  forward  Mach 
cone  at  any  point.  For  these  problems,  Woodward  (Reference  9)  recently  suggested  a 
third  order  singularity  which  would  still  allow  the  panel  methods  to  be  utilized. 

From  the  above  discussion,  we  could  see  that  the  present  method  of  integral  equation 
could  handle  fairly  complicated  flows  within  the  limitations  of  moderately  non-linear 
flows,  particularly  for  off-body  points  velocity  computations. 
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WEAK  VISCOUS  -  INVISCID  INTERACTION 

3.1  Viscous  Flow  Modelling  in  Panel  Methods 

The  simplest  and  most  popular  method  of  viscous  correction  to  an  inviscid  flow 
computation  is  discussed  in  this  section.  In  this  process  the  matching  line  is  not  used 
directly.  Instead,  a  displacement  surface,  deduced  from  the  integral  mass  conservation 
property  of  the  boundary  layer  replaces  the  real  surface  for  the  inviscid  part  of  the 
calculation.  This  displacement  surface  is  treated  as  a  streamline  of  the  inviscid  flow. 
The  implicit  assumptions  in  this  process  are: 

(a)  The  displaced  surface  is  a  real  streamline. 

(b)  The  pressure  computed  at  the  displacment  surface  in  the  inviscid  calculation  is  the 
same  as  the  pressure  which  drives  the  boundary  layer. 

These  assumptions  are  valid  as  long  as  the  viscous  effects  are  weak.  To  link  viscous  and 
inviscid  influences  through  displacement  thickness,  there  are  two  procedures: 

(«)  Actually  displacing  the  body  surface  by  the  addition  of  the  displacement  thickness 
locally. 

(ii)  Implicitly  represent  the  displacement  effect  by  the  surface  blowing  or  transpiration 
method,  by  imposing  artiflcal  boundary  conditions  on  that  body  which  simulates 
the  displacement  surface  as  a  streamline  of  the  inviscid  flow. 

Method  (ii)  is  recommended  whenever  iterations  are  involved. 

3.1.1  Displacement  Thickness  Method 

The  idea  of  adding  the  displacement  thickness  to  the  inviscid  flow  to  incorporate 
the  weak  viscous  effects  could  be  demonstrated  as  follows  (  Reference  10). 

The  first  order  outer  expansion  in  the  asymptotic  expansion  of  the  equations  of 
motion  with  the  viscosity  as  a  small  parameter  is  : 
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(11) 
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dz2  dy2 

where  rh\  is  stream  function  and  Bfl  is  related  to  the  vorticity  . 


Since  the  vorticity  is  a  function  of  ipi  only,  it  must  be  constant  along  streamlines. 
For  the  problem  of  flow  over  a  flat  plate  with  a  uniform  upstream  irrotational  flow  the 
solution  becomes. 

i)\{x,y)  —  Uy  (12) 

which  means  a  uniform  parallel  stream  with  a  velocity  U. 

The  first  order  inner  solution  for  the  semi-infinite  flat  plate  with  Y  =  >/ Re y  .  Re 
being  the  Reynold’s  number,  is  obtained  from  the  equation  of  motion  as: 

Q  ^ 

^lrrr V'l,,— )^ir  =0  (13) 

with  the  associated  boundary  conditions  , 


^i(ijO)  =  0  pir(x,0)  =  0 


and 


u-’iy(:r.  oo )  =  U 


(14) 


where  the  subscripts  x  and  Y  again  mean  derivatives  with  respect  to  the  same 
variables. 


Using  one  parameter  Group  Transformation,  we  can  obtain  the  invariants  of  the 
transformation  as  : 

V  =  and  /i(r?)  —  ~z  (15) 

s/2z  V2x 

and  substituting  these  transformations  into  the  differential  equation  (13)  gives  the  clas¬ 
sical  Blassius  boundary  layer  equation  , 


fT  +  Ml  =  0  (16) 

with  the  boundary  conditions  , 


/i(0)  =  /\(0)  =  0  'ind  ?,{<*)  =  U. 


(IT) 


and  the  primes  denote  differentiation  with  respect  to  rj 


The  second  order  outer  expansion  for  the  above  problem  is  governed  by 


d2^ 2  ,  92<u 2  _ 
dx 2  dy2 

The  boundary  conditions  near  the  plate  are  obtained  by  the  matching  condition, 


(18) 


^2(2.0)  =  0  for  x  <  0 


V>2(z,0)  =  —  Ltx^ociXibxy  —  ipi)  =  —  \flx3  for  x  >  0 
The  condition  far  upstream  is 


(19) 


02  {x,y)  =  0(y)  (20) 

This  is  the  linearized  problem  for  potential  flow  past  the  displacement  parabola  Y  = 
,  where  0  =1.21678  .  The  second  order  inner  expansion  could  be  obtained  as 
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dY 


Vlyyy 
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=  0 


(21) 


and  the  associated  boundary  conditions  , 


ifj2{x,  0)  =  02y(z-  0)  =  0,  and 


V2y(-E,y)  =  0  as  y  h  oc  (22) 

Hence  ,  the  complete  second  order  approximation  for  the  flow  past  a  semi-infinite  flat 
plate  is 

outer  solution  : 


%>(x,Y  ;Re)  =  UY - —dxlJi{x  +  iY) 

VRe 


(23) 
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inner  solution: 


xp{x,Y;Re) 


VRey 


+0(s» 


(24) 


From  the  above  explanations,  we  can  see  that  the  method  of  adding  the  displace¬ 
ment  thickness  by  any  of  the  two  methods  to  account  for  the  viscous  effects  in  the  case 
of  wing-body  configuration  is  a  good  first  order  approximation  as  long  as  the  viscous 
effects  are  weak.  Further,  when  the  process  is  continued  once  more,  i.e.,  obtaining  the 
second  order  inner  solution  and  deriving  the  third  order  outer  solution,  the  contribution 
of  this  order  approximation  to  the  final  solution  is  insignificant.  Hence,  we  can  stop  the 
viscous/inviscid  interaction  calculation  after  one  iteration  to  arrive  at  a  good  approxima¬ 
tion  (Reference  11  and  12  ).  Also,  Reference  11  shows  that  a  Strip  Theory  approximation 
by  solving  the  two-dimensional  compressible  boundary  layer  equations  is  adequate  for 
the  inclusion  of  viscous  effects  compared  to  a  three-dimensional  compressible  boundary 
layer  equations  with  small  cross  flow  assumption. 


3.1.2  Surface  Blowing  Method. 

Since  the  displaced  surface  is  assumed  a  streamline,  the  boundary  condition  on  6* 
is  that  in  inviscid  calculations  the  flow  direction  is  tangential  to  the  displaced  surfaces. 
With  subscripts  5*  and  s  meaning  the  quantities  on  the  displaced  and  on  the  original 
surf  aces,  respectively, 


(25) 


t >$•  _ 

us*  ds 

To  apply  transpiration  method,  we  can  write  is *  as  a  series  expansion  of  the  velocity  on 
the  surface  with  5*  as  a  small  parameter  and  with  the  usual  definitions  of  u  and  v  being 
the  compressible  velocities  in  the  s  and  n  directions,  along  and  normal  to  the  surface, 
respectively,  i.e., 


But  from  continuity, 


+  5*— |» +0(5  ) 


(26) 
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using  (27)  and  (25)  in  (26) 


08 *  *<9u 

-  u‘-g7  + s  si  '* 


Again,  writing  as  a  series  expansion, 


,  c*0u  .  ,  ,  *2. 

=  +  ^  U+0(^  ) 

On 


and  using  irrotationality  criterion, 


and  substituting  it  in  (29) 


On  _  dv 

On  Os 


.  ,  c*0v  .  . 05  c*0u 

vs  =  (us  +  8  5  U  (3°) 

This  is  the  surface  blowing  normal  velocity  in  terms  of  the  velocities  on  the  real  surface 
which  can  be  calculated  from  the  same  influence  coefficient  matrix  for  each  iteration.  The 
gradients  can  be  calculated  using  central  differencing  with  the  5*  values  being  computed 
using  a  finite  difference  or  integral  boundary  layer  calculation  (Reference  13)  procedure. 
As  the  trailing  edge  is  approached.  is  an  order  of  magnitude  greater  than  us  and 
hence  5*^  |s  is  not  negligible. 

The  pressure  coefficient  cp  is  given  by 

cp  =  1  (<78  )g* 

=  1  -  (uf.  +  v%.) 

Cp  =  1  -  (us  +  6*J-  |s)2  -  (tt,  -  d'jj-  |s)2  (31) 

The  present  scheme  for  viscous-inviscid  interaction  procedure  with  equivalent  blow¬ 
ing  approximation  uses  a  three-point  second  order  finite  difference  parabolic  fit  for  the 
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blowing  velocity 


Defining,  q  —  ue6*, 


v *  -  K”*4') 


(32) 


( ^1\  —  I*  ~  Sbzzi 
\ds  )  k  Sk  —  Sk~  i 


1  — 


Sk  —  <S*  —  1 
Sfc  +  1  —  -S*  —  ! 


+ 


qk+ 1  —  ?fc 


5fc-(-l  Sfc 


1  - 


•S  fc  — l  —  Sk 
S*+l  ~  Sk —  I  J 


for 


*  =  1,2, ---.(.V—  1) 


and 


f  dq\  _  (s,y  —  ajy+i)<7iV— 2  _ sy  — 

\ds  )  N  (  5  ,V — 1  —  S.v— 2)(s,V  —  Sy  — 2)  (S,V—  1  —  S.V— 


—  Sy_ 2 


>)(SjV  —  S  ;V - 1  ) 


2  Sy  —  Sy— 1  —  Sy— 2 

(■s.v  —  s.v— 2)(s,v  —  -Sy— 1) 


?iV 


(33) 


?V— 1 


(34) 


for  the  last  panel. 

3.1.3  Kutta  Condition 


Inviscid  Flow 
Explicitly: 

(a)  Equate  the  upper  and  lower  surface  velocities  at  or  near  the  trailing  edge. 

(b)  Specify  a  velocity  component  usually  (but  not  necessarily  zero)  normal  to  the 
extended  camber  in  the  wake. 

Implicity: 

(a)  The  basic  formulation  of  the  mathematical  model  inherently  includes  a  condition 
which  yields  a  smooth  flow  at  the  trailing  edge. 

Viscous  flow 


The  rule  of  stagnation  pressures  (obtained  from  potential  flow)  at  the  non-cusped 
trailing  edge  to  control  circulation  is  not  correct.  But  from  boundary  layer  assumption 


and  the  matching  line  assumption,  we  can  imply  that  the  pressure  outside  the  displace¬ 
ment  surface  is  the  same  as  that  on  the  body.  Hence,  the  condition  to  control  circulation 
is  that  pressures  are  equal  on  the  displacement  surface  on  both  sides  at  the  trailing  edge. 
This  criterion  transforms  in  the  case  of  the  surface  blowing  case  as 

%2:  +  vbi  =  (35) 

(i;  -  vh) 

-  us;  =  7 — — - r  (36 ) 

which  is  similar  to  the  inviscid  case  of  equal  velocities  when  Us *  =  UgJ-  Equation  (36! 
can  be  written  as: 


=  u,;  +  [v„  -I'lj;  l«.)!  - 


(v-su  +  K  ff  |su)  +  («s,  -h  Uj 


c*  dv  .  c*du 

(3.) 

Note  that  all  the  terms  containing  6*  are  on  the  RHS  and  consequently  the  influence 
coefficient  matrix  remains  unaltered.  The  subscripts  u  and  /  for  5*  and  •?  denote  the 
upper  and  lower  surfaces  for  the  last  panel  near  the  trailing  edge.  No  special  treatment 
has  been  made  for  the  Kutta  condition  except  that  the  usual  condition  of  equal  pressures 
on  the  centroid  of  the  last  panel  on  the  top  and  bottom  displaced  surfaces  is  imposed. 

The  boundary  layer  equations  are  solved  by  the  Keller's  Box  Method  (Reference  13 
and  14).  The  flow  quantities  are  solved  at  the  center  of  each  mesh  with  the  pressure 
gradient  obtained  from  the  external  inviscid  velocities.  We  note  that  the  displacement 
thickness  is  the  only  boundary  layer  parameter  that  affects  the  modified  potential  flow. 
Special  care  has  been  taken  to  determine  the  stagnation  point  of  each  lifting  section. 
Computationally,  this  location  is  found  by  taking  the  dot  product  of  the  velocity  at  the 
mid-points  of  each  panel  with  the  average  of  the  unit  tangent  vectors  of  the  two  N-lines 
bordering  the  panel.  This  dot  product  changes  sign  exactly  at  one  point  in  each  lifting 


strip.  The  desired  stagnation  point  is  obtained  by  interpolation  with  respect  to  values  of 
the  dot  product  between  the  mid-points  on  either  side  of  the  sign  change.  The  boundary 
layer  calculation  can  be  carried  out  for  both  laminar  and  turbulent  cases.  The  turbulent 
boundary  layer  calculations  are  done  using  a  mixing  length  model.  The  calculations 
predict  the  separation  point  .  The  present  code  can  automatically  incorporate  any  of  the 
above  mentioned  two  viscous  effects  at  only  one  angle  of  attack  at  a  time,  compared  to 
the  possibility  of  calculating  the  potential  flow  for  up  to  ten  angles  of  attack. 

For  small  angles  of  attack  and  thin  wings  the  boundary  layer  effects  are  small. 
A  further  modification  of  the  code  to  incoporate  the  boundary  layer  computation  on 
axi-symmetric  body  is  effected.  This  requires  the  non-lifting  body  panel  M-lines  and 
N-lines  to  be  interchanged  in  direction  compared  to  the  potential  part.  The  panels  are 
redistributed  from  tip  to  nose  on  the  bottom’ surface  followed  by  nose  to  tip  on  the 
top  surface.  The  necessary  modifications  in  the  panel  surface  normal  velocity  vector 
calculations  and  the  velocity  vector  computations  are  made.  No  special  treatment  for 
stagnation  point  determination  has  been  introduced,  since  most  non-lifting  bodies  under 
study  have  a  fairly  pointed  nose  shape.  Another  justification  for  this  is  the  fact  that  a 
slight  deviation  of  the  stagnation  point  location  is  not  expected  to  alter  the  store  flow 
field  appreciably.  A  Fortran  listing  of  all  these  modifications  is  included  in  Appendix  C. 

3.1.4  WAKE  EFFECTS 

The  presence  of  the  wake  has  two  distinct  contributions  on  the  wing.  They  are 
essentially 

(i)  Displacement  effect. 

(ii)  Curvature  or  pressure  discontinuity  effect. 

The  displacement  effect  of  the  wake  is  precisely  analogous  to  that  of  the  boundary  layer 
over  the  solid  surface  of  the  wing.  Thus  it  can  be  represented  either  by 

(a)  source  distribution  over  an  infinitesimally  thin  surface,  leading  to  a  discontinuity  in 

velocity  normal  to  the  surface  given  by, 
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(b)  as  a  streamtube  of  thickness  5*w. 
where  the  subscript  w  refers  to  the  conditions  on  the  wake  surface. 

In  either  case  the  normal  first-order  definition  of  8*w  will  suffice,  since  the  effect  of  the 
wake  on  the  pressure  distribution  of  the  wing  is  small. 

To  compute  &*w  on  the  wake  using  the  panel  method,  the  following  procedure  might 
be  used: 

(1)  During  the  inviscid  computation  process,  distribute  off-body  points  close  to  the 
approximate  wake  surface  and  obtain  the  inviscid  velocites  at  those  points. 

(2)  In  the  boundary  layer  calculations,  modify  the  airfoil  shape  locally  to  include  this 
wake  surface  and  hence  obtain  the  5*w. 


For  an  inviscid  flow,  and  approximately  in  viscous  flow  (Figure  2), 


where  k  is  curvature  of  the  streamline  and  p  is  the  density.  Integrating  across  the  wake, 
we  have 

Pn  —  Vi  =  ~  /  kpu2dn 


Ph  —  Pi  =  —  J  kpu2dn 

Pu  —  Pi  =  — kpwu2Sw  +  tf* } 


where  5W  and  9W  are  the  boundary  layer  and  momentum  thicknesses, respectively.  Again, 
considering  the  pressure  jump  across  the  wake  in  the  sense  of  source  distribution, 


Ap  =  {pu  —  Pi)  +  kdwu2Jw 
Ap  =  kpu2m(i 5*  + 


since  A p  is  small,  it  could  be  written  as 


iu  =  + 


The  effect  described  above  is  similar  to  that  of  a  Jet  flap,  but  of  opposite  sign,  with  the 
jet  momentum  coefficient  replaced  by  {<5^  +  0W). 


The  determination  of  the  exact  wake  shape  is  fairly  complex  and  hence  within 
the  present  order  of  approximations,  we  can  make  certain  simplifying  assumptions.  An 
effective  way  to  determine  the  shape  of  the  wake  is  to  require  that  the  wake  surface 
originates  from  the  trailing  edge  and  that  its  tangent  bisects  the  total  velocity  vectors 
on  the  upper  and  lower  side  of  the  wake  surface  i.e., 


f  vnm  A  _  ( 

V  U8W  J  u  V  USmJ( 


(42) 


In  agreement  with  the  boundary  layer  theory,  for  the  determination  of  the  wake  shape 
alone  we  impose  zero  pressure  jump  across  the  wake  (curvature  effect  is  negligible),  which 
implies 


(uOu  +  (“sJu  =  (t'njf  +  (Ue„)u 


(43) 


from  which  follows  that 


=  (VnJe  and  («».)»  =  (“>.)<  on  the  wake  (44) 

Thus,  the  tangential  velocity  is  continuous  across  the  wake  and  the  normal  velocity  is 
symmetric,  and  thus  discontinuous,  with  respect  to  the  wake.  Thus,  the  wake  effect 
could  be  represented  by  a  source  distribution  on  the  wake  only. 

For  non-lifting  bodies,  the  effect  of  wake  is  necessary  to  be  included.  The  wake  of 
the  non-lifting  blunt  based  body  is  quite  complex  because  of  the  base  flow.  But  within 
the  order  of  approximations  incorporated  in  the  other  parts,  the  wake  of  a  non-lifting 
axisymmetric  body  could  be  represented  as  an  effective  body  of  equivalent  base  shape 
for  inviscid  calculations.  The  shape  of  the  equivalent  base  could  be  approximated  as 
a  backward  facing  step  in  two  dimensions.  Previous  experiments  conducted  in  UTSI 
(Reference  15)  have  yielded  that  inviscid  streamline  for  the  backward  facing  step  reat¬ 
tached  the  sting  behind  the  base  of  the  step  at.  approximately  six  base  heights  from  the 
base.  This  relation  was  utilized  when  repanelling  the  non-lifting  bodies  for  the  inclusion 
of  viscous  effects. 


3.1.5  Separation  Effects 


Inclusion  of  weak  separation  effects  in  an  overall  sense  is  another  step  toward 
improving  an  airfoil  or  wing  viscous-inviscid  interaction  computational  scheme.  The 
separation  effects  could  be  due  to  strong  or  weak  viscous  inviscid  interactions.  The 
flow  separation  due  to  strong  interactions  like  shock  boundary  layer  interaction  is  very 
complex  to  incorporate.  The  weak  interaction  induced  separation  could  be  successfully 
modeled  in  the  following  way. 

The  near  flow  field  of  the  airfoil  is  divided  into  four  different  regimes  (Figure  3). 
The  conventional  inviscid  panel  methods  could  be  applied  at  Region  1.  For  Region  2. 
the  integral  or  differential  boundary  layeT  equation  approximations  are  adapted.  For  the 
free  shear  layers,  Region  3,  vortex  sheets  of  constant  density  along  their  length  starting 
from  the  upper  and  lower  separation  points  (the  lower  surface  separation  point  being 
taken  as  the  trailing  edge  for  positive  angle  of  attack).  The  shape  of  the  vortex  sheet 
is  not  known  a  priori,  and  it  is  iteratively  determined  as  a  streamline  of  the  resultant 
flow.  The  vortex  strength  on  the  upper  and  lower  surface  sheets  are  set  to  be  equal. 
Thus,  the  modified  airfoil  shape  is  determined  for  the  next  inviscid  computations.  In 
Region  4,  the  reduced  total  head  must  be  taken  into  account  when  calculating  the  body 
pressures  at  points  downstream  of  separation.  The  two  basic  assumptions  embedded  in 
this  categorization  of  the  regions  are: 

(1)  The  boundary  layer  and  the  free  shear  layers  are  thin  and  hence  could  be  represented 
as  slip  surfaces  which  are  streamlines  across  which  there  exists  a  jump  in  velocity. 

(2)  The  wake  does  not  have  significant  vorticity  and  has  constant  total  pressure  (lower 
than  the  free  stream  total  pressure)  and  hence  a  potential  flow  region. 

The  usual  boundary  conditions  of  V.n  =  0  for  potential  flow  and 


for  viscous-inviscid  interaction  are  applied. 
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The  following  approximate  relations  are  empirically  determined  by  earlier  investigators 
(Reference  16)  and  could  be  judiciously  used  (see  Figure  3  for  symbols) 


0.192  ( 0.72  Y2  , 

WL  =  W„  +  1.118  +  (40) 

For  the  upper  shear  layer, 

Vm  =  ~(V0Uter  +  Vinner)  =  average  velocity  in  the  layer  (46) 

where, 

V0uter  —  Vm  “f“  ;^7u 

Vinner  —  Vm  “  7u  a^d 

T u  ==  V0uter  Tjnner 

is  the  vorticity  in  the  upper  shear  layer  (and  vice  versa  for  the  lower).  The  zero  static 
pressure  drop  across  the  shear  layer  implies  that  the  total  pressure  jump  across  the  shear 
layer  could  be  calculated  as 

(47) 

with 

AH  —  Hinner  Hotter 


AH  —  Pinner  “1“  2^”^  Pouter  ~t~  ~P^'m  "f" 


Since  the  wake  is  assumed  to  have  constant  total  pressure,  the  jump  in  total  pressure 
across  the  free  shear  layer  is  the  same  everywhere  in  the  wake. 

The  modified  pressure  coefficient  is  given  by 


—ta 


+ 


AH 

Qoc 


(48) 


where  Vx  and  qoo  are  the  free  stream  velocity  and  dynamic  pressures,  respectively. 
The  above  discussion  is  appropriate  only  for  weak  separation  which  occurs  at  the  trailing 
edge  of  wing  at  moderate  angles  of  attack.  Strong  separation  induced  by  shock  boundary 
layer  interactions  and  by  high  angles  of  attack  are  more  complex  (Reference  17). 


SECTION  rv 


INTERNAL  FLOW  COMPUTATION 

The  results  obtained  with  the  method  of  flow-field  velocity  computation,  as  dis¬ 
cussed  in  Section  II, were  in  reasonable  agreement  with  the  experimental  data  for  most 
configurations  excepting  those  which  included  field  points  in  the  gap  region  that  appears 
between  the  wing  and  the  store  (Configurations  101  and  102).  These  configurations  es¬ 
sentially  simulated  stores  shortly  after  separation.  The  transonic  flow  in  these  regions 
should  contain  shocks,  as  observed  in  experimental  data,  and  strong  viscous  inviscid 
interactions.  The  shocks  in  these  regions  were  smeared  off  in  the  linear  theory,  and  the 
non-linear  equation  method  predicted  wrong  location  and  strength  of  the  shocks.  This 
region  should  have  to  be  analyzed  by  a  much  more  elaborate  method  with  the  represen¬ 
tation  of  the  physics  of  this  complex  flow  through  an  adequate  mathematical  modelling. 
Hence,  in  the  hierarchy  of  the  approximate  analysis  of  the  flow  field,  this  gap  region 
has  been  investigated  independently  by  a  finite  difference  method  for  solving  the  two- 
dimensional  full  Navier-Stokes  equations  (the  two-dimensional  section  being  the  vertical 
plane  passing  through  the  store  centerline  where  the  experimental  data  was  available). 
The  Navier-Stokes  equations  were  written  in  characteristics  form.  MacCormack's  ex¬ 
plicit,  time  marching,  factored,  predictor-corrector,  finite  difference  scheme  (Reference 
18)  has  been  utilized  to  compute  the  internal  flow  between  the  two  walls  (underneath 
the  wing  and  above  the  store).  Because  of  the  inherent  resistance  of  the  streamline  to 
change  their  curvature  at  transonic  flow,  the  two-dimensional  approximation  is  relevant. 
The  approximate  inlet  and  exit  boundary  conditions  at  the  appropriate  mesh  points  are 
obtained  from  the  external  flow  computations,  treating  those  points  as  off-body  points. 

The  two-dimensional  Navier  Stoke's  equations  are  written  as: 


0 u  .  0u  .  du 

m+'lr,  +  vTy  + 


1  dp  1  o  L  .  „  .du  .  s  fli/1  ,  1  0  r  ,6v  .  du. 1 
-poi  “  -pei r  +  +  X^J +  pdA^di  +  ^}J 

0v  .  dv  ,  Dv  ,  i  dp  1  d  L.  ,  .  ,0t>  ..fltil  1  d  \  0v  flu] 
~0t  ^  U0j  **"  V0y  pdy  “  pdy[~  ^0y^~  0*]  pdx{  Bx  fly] 
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where  p  is  the  density,  p  is  the  pressure,  T  is  the  temperature,  u  and  v  are  the  velocity 
components  in  x  and  y  directions,  a  is  the  speed  of  sound,  R  is  the  gas  constant,  p  and 
X  are  the  first  and  second  coefficients  of  viscosity,  respectively;  ^  is  the  ratio  of  specific 
heats,  k  is  the  thermal  conductivity,  x  and  y  are  the  space  coordinates,  t  is  the  time. 


The  equations  are  solved  using  a  predictor-corrector  explicit  MacCormack's  scheme. 
The  predictor  step  uses  a  backward  differencing  and  the  corrector  step  uses  a  forward 
differencing  to  the  convection  terms.  Central  differencing  is  used  for  the  diffusion  terms 
in  both  predictor  and  corrector  steps.  The  scheme  is  stable  provided  that  the  time  step 
increment  satisfies  the  stability  criterion, 


S  +  +  ay/dx*  +  3^ 


where  A  is  a  parameter  like  Courant  Fredrick  Levy  (  CFL  )  number  which  neeeds  to 
be  set  less  than  1  for  the  cases  of  supersonic  calculations  with  shock.  A  fourth  order 
artificial  viscosity  and  second  order  artificial  thermal  conductivity  terms  are  added  to 
dampen  the  numerical  oscillations  near  shocks.  The  time  split  procedure  allows  the 
marching  in  two  different  directions  independently  by  introducing  a  term  which  is  of 
the  same  order  as  the  truncation  error.  The  steady  state  solutions  are  obtained  by  the 
asymptotic  unsteady  calculations. 

The  solution  procedure  automatically  computes  the  initial  two-dimensional  flow- 
surface  by  assuming  one-dimensional  approximation  based  on  the  pressure,  velocity,  and 
density  values  specified  at  the  inlet  of  the  convergent  divergent  section.  The  boundarv 
conditions  of  impermeable  and  either  slip  or  no-slip  walls,  together  with  the  exit  pressure 
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conditions,  are  satisfied  to  obtain  the  final  solution  surface.  The  truly  external  flow 
calculation  would  yield  maximum  velocities  at  those  regions  which  are  closest,  i.e.,  at 
the  minimum  area  regions,  whereas,  the  internal  flow  calculation  would  have  the  throat 
(  minimum  area  )  choaked.  The  pressure  distribution  upstream  readjusts  itself  so  as  to 
yield  sonic  velocity  near  the  throat.  There  may  be  a  further  acceleration  to  supersonic 
speed  downstream  and  finally  recover  the  exit  pressure  through  a  near  normal  shock. 
For  the  present  problem  involving  the  flow  in  the  gap  between  the  pylon  or  wing  and 
the  store  is  not  essentially  two-dimensional  or  completely  internal  flow. 

The  basic  YNAP  (Viscous  Nozzle  Analysis  Program)  scheme  described  above  does 
pose  difficulties  to  capture  shock  at  transonic  speeds.  The  scheme  was  essentially 
developed  for  fully  subsonic  or  fully  supersonic  and  plug  nozzle  flows.  The  capturing 
of  shocks  should  be  done  by  carefully  adjusting  the  artificial  viscosity.  Originally, 
the  computed  flow  either  continued  to  accelerate  to  supersonic  speeds  or  smoothly 
decelerated  subsonically  downstream  of  the  sonic  points  for  the  cases  of  supersonic  or 
subsonic  exit  boundary  conditions,  respectively.  In  order  to  modify  the  scheme  to  obtain 
a  solution  with  shock,  an  indirect  approach  to  fit  the  shock  by  a  forward  and  backward 
sweep  process  has  been  adapted.  During  the  forward  march  normal  shock  was  assumed 
to  be  present  at  all  supersonic  stations  and  the  Rankin-Hugoniot  jump  relations  have 
been  imposed  at  all  those  stations.  The  downstream  Mach  numbers  across  the  normal 
shock,  at  each  supersonic  station  is  then  computed.  A  reverse  sweep  from  the  exit 
station  with  the  approximate  switched  inlet  boundary  conditions  is  then  made  with  the 
modified  wall  and  centerbody  shapes.  These  modified  shapes  are  obtained  by  adding 
the  displacement  thicknesses  (available  from  the  previous  forward  sweep)  to  the  original 
shapes.  The  points  where  the  matches  with  the  Mach  number  obtained  during  the 
backward  sweep  M~  determines  the  shock  location.  Figures  4,  5.  and  6  show  the  above 
idea  pictorially.  The  shock  locations  determined  from  the  procedure  discussed  above 
compare  well  with  those  observed  in  the  experimental  results  for  both  the  configurations. 
Figure  7  and  S  show  the  direct  computation  of  the  shock  location  for  the  above  two  cases 
by  carefully  adjusting  the  artificial  viscosity  terms.  The  same  values  of  the  artificial 
viscosity  terms  are  used  for  the  other  configuration  internal  flow  computations  as  well. 
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SECTION’  V 


COMPUTATION  OF  STORE  PRESSURE  DISTRIBUTION 

The  store  pressure  distribution  computation  procedure  described  herein  incorporates 
the  attractive  features  of  the  different  methods  of  solution  of  the  transonic  multibody 
problem.  The  flexibility  of  the  panel  code  for  complex  geometries  and  the  accuracies 
of  the  finite  difference  codes  are  fused  together  to  yield  an  approximate  method  which 
does  not  involve  very  elaborate  computations.  The  formulation  of  the  problem  and  the 
solution  procedure  for  store  pressure  calculation  by  the  RAXBOD  (Relaxation  solution 
of  AXi-symmetric  BODv)  finite  difference  scheme  is  outlined  below  (Reference  19). 


The  exact  equation  for  the  disturbance  potential,  <p,  for  axi-symmetric  compressible 
flow  can  be  written  in  orthogonal  curvilinear  coordinates  £,  7  as: 
f  u2\  1  d  (  1  d®\  0  uv  d2<t>  f  u2\d20  ["  t  j (  cos# 

4wt  +  — 1^  =  0  08, 

Ha 2  r  H  dE, 


where  6  and  k  are  the  angle  (measured  counterclockwise  from  the  axis  of  symmetry)  and 
curvature  of  the  reference  coordinate  surface, respectively,  H  =  1  -j-  krj  is  the  metric, 
a  is  the  local  speed  of  sound,  u  and  V  are  the  £,  7  components  of  the  total  velocity,  q, 
and  r  is  the  radius  from  the  axis. 

Bo 

—  H(u  cos#) 
d(t>  ,  .  „ 


=  v  -f  sin( 


The  coordinate  system  is  built  in  a  special  way  so  as  to  have  7  =  0  on  the  body  surface 
throughout  and  to  remove  the  difficulties  associated  with  the  fully  body  fitted  coordinate 
system  in  the  concave  part  of  the  body  surfaces.  Hence,  in  the  forebody  region,  the  body 
normal  coordinates  are  used  upto  the  first  horizontal  tangent  and  beyond  that  point  a 
sheared  cylindrical  system  has  been  introduced. 


A  convenient  normal  coordinate  stretching  function  has  been  applied.  It  is 
represented  as 


Cd|  Qj 


(49) 


Ay 

n  =  - 

(i  -  y)a 

where  rj  is  the  physical  coordinate  normal  to  the  body  and  y  is  the  computational 
coordinate  which  varies  from  zero  at  the  body  surface  to  one  at  infinity.  The  constant 
A  could  be  used  to  control  the  physical  step  size  at  the  body.  A  =  (^)y=o  and  for 
a  given  value  of  .4,  the  exponent  a  is  used  to  control  the  size  of  the  last  finite  value 
of  tj  (the  larger  a  is,  the  farther  the  points  are  away  from  the  body).  The  tangential 
coordinate  stretching  is  based  on  the  physical  arc  length  along  the  reference  surface. 


5  =  +  U  -  i)|.4  +  B(i  -  i)2 

(50) 

where  A  and  B  are 

determined  by  specifying  (37^  and  requiring  that  £ 

at  X  =  1.  Hence. 

3£max  ~  (§)xi=0 

(51) 

A  =  - 

2 

and 

B  —  4  [£rnas  —  A) 

■  (52) 

For  open  bodies  the  tangential  coordinate  stretching  is  divided  into  two  regions  with  the 
physical  location  of  the  dividing  point,  xm.  The  stretching  Junction  used  for  the  region 
from  the  nose  up  to  xm  is 


£  =  \,x  -j-  X2 -)-  X3J5  -+•  X4Z'  ,  0  <x<xm  (53) 

and  in  the  region  from  xm  to  unity,  the  stretching  function 

c  _  r  \  j,  ^  Xrn  1  rr.  )  ^  ~  ^  1  /  r  .1  1 

s  —  “i  "  ,  Xm  _~.X  <  l  (o4) 

(1  —  x) 

is  applied.  The  coefficients  in  these  expressions  are  determined  by  specifying 
£mt(3j)x=o,  and  jj|  be  continuous  at  x  =  xm. 

In  the  sheared  cylindrical  coordinates,  the  potential  equation  is 

f  u2\<920  Jdrb[]  u2\  uv  d2® 

V  a2)d£2  ~[  dx  V  a2/  a2  d^drj 
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+  +  :  F  =  0  (55) 

ox i  a z  r\orj 

where  £  is  again  identified  with  the  axial  coordinate  x,  and  r]  is  a  transformed  radial 
coordinate  such  that  r]  =  0  on  the  body  surface,  and  rj,  is  the  local  body  radius.  In  the 
sheared  cylindrical  system, 

do  drb  do 

—  =  u  -f  — —  v  —l,  and  —  =  v  (o6) 

d£  dx  dr] 

The  boundary  conditions  are 


T]-*0 0 


r?  =  0 


v  —  0, 


— —  =  sin  (9  for  x<xm 

drj 


auu 

”  r,  =  ^Twh+Ti)  r°'  ^  m 

The  computational  scheme  incorporates  a  rotated  differencing  at  supersonic  points, 
where  u2,v2  <  a2  <  u2  -f-  v2 ,  to  avoid  the  possible  difficulties  that  would  arise 
during  the  solution  process  because  of  the  misalignment  of  the  sonic  line  with  the  r] 
coordinate  surface.  The  rotated  differencing  scheme  replaces  the  highest  derivative  part 
of  the  differential  equation  (48)  from  £,  T]  coordinates  to  5,  n  coordinates  which  are  nearly 
along  and  normal  to  the  velocity  vector  direction  locally  by 


0ss  “I"  <^nn  =0 
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where 


1  I" v2  d  f  1  dC>\  uu  d24>  2<920 

'ss  ~  q2  [i7  Wi;\H  h)  +  h  dtdr]  +  v  <V 

_  1  v2  d  {  1  ^uvdct) 
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uv  d<t>  od20 


for  the  bodv-normal  coordinates  and 
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Further  details  of  the  above  procedure  could  be  obtained  from  Reference  20.  The  code 
incorporates  a  successive  mesh  refinement  from  a  coarse  to  a  fine  mesh  as  many  sequential 
steps  as  desired  within  the  purview  of  storage  space  and  CPU  time  limitations. 


The  present  store  pressure  computation  procedure  follows  the  steps  mentioned 


below: 


(1)  Compute  the  store  only  and  the  entire  configuration  on-store  pressure  distributions 
using  the  panel  method. 

(2)  Calculate  the  interference  effects  by  subtracting  the  store  alone  on-body  pressures 
locally  from  that  of  the  entire  configuration  store  on-body  pressures. 

(3)  Compute  the  store  only  on-store  pressure  distribution  using  the  R.VXBOD  finite 
difference  code.  Note  here  that  a  linear  or  higher  order  interpolation  in  velocity  is 
needed  to  find  the  pressures  at  the  same  points  as  those  of  the  panel  method. 

(4)  Add  the  pressure  coefficient  due  to  interference  and  the  store  alone  on-body  finite 
difference  computations  to  yield  a  corrected  on-store  pressure. 


The  above  procedure  is  questionable  in  the  sense  that  the  linear  combination  of  the 
flow  quantities  is  attempted  in  a  non-linear  flow  (both  due  to  geometry  and  due  to  the 
equations  of  fluid  motion).  A  more  reasonable  approach  would  be  to  build  in  correlations 
between  the  store  alone  and  wing-body-store  combination  on  store  pressures  and  between 
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the  panel  and  finite  difference  methods  at  the  same  transonic  Mach  numbers  for  the 
store  alone  case.  Several  such  efforts  have  been  tried,  and  at  the  time  of  reporting  a 
comprehensive  correlation  was  not  arrived  at.  In  retrospect,  the  results  of  the  procedure 
discussed  formerly  have  only  been  presented.  It  could  be  seen  that  the  above  scheme 
did  not  improve  panel  method  prediciton  very  much.  A  much  more  elaborate  way 
of  incorporating  the  finite  volume  codes  (FLO-27)  for  wing-body  configurations  and 
panel  method  for  the  entire  flow  field  is  being  investigated,  and  it  is  speculated  to  yield 
improved  predictions. 
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SECTION  VI 


RESULTS  AND  DISCUSSION 

In  this  section  the  resuits  of  the  calculation  procedures  discussed  in  the  previous 
sections  for  different  configurations  are  presented.  The  present  study  has  been  mainly 
directed  towards  improving  the  results  obtained  in  the  previous  approaches  and  arriving 
at  a  unified  and  simple  procedure  for  the  analysis  of  complex  wing-body-pylon-store 
configurations  at  transonic  speeds. 

The  experimental  data  base  for  the  present  theoretical  comparison  comprises  the 
data  obtained  from  the  Arnold  Engineering  Development  Center  (AEDC)  lfoot  x 
Ifoot  transonic  wind  tunnel  (IT)  Laser  Velocimeter  flow-field  survey  experiments  for 
configurations  (a)  shown  in  Appendix  A,  and  from  the  AEDC  lfoot  x  lfoot  tran¬ 
sonic  wind  tunnel  (4T)  5  hole  pressure  probe  flow-field  survey  for  Configuration  24  of 
configurations  (b)  given  in  Appendix  A.  The  configurations  (a)  have  wall-mounted  45- 
degree  swept  wing  with  Aspect  Ratio  (AR)  of  2.73  having  NACA  0005-31  airfoil  sections 
mounted  with  two  different  store  shapes.  The  free  stream  Mach  numbers  of  test  were 
0.92,  0.975,  and  1.025  with  wing  angles  of  attack  of  0  degree  and  2  degrees  at  Reynolds 
number  of  about  3  million  per  foot.  Configuration  001  is  similar  to  Configuration  221. 
but  with  the  inclusion  of  a  fuselage-mounted  wing  instead  of  the  wall-mounted  wing. 
Configuration  001  has  been  tested  at  a  subcritical  Mach  number  0.7.  Configurations 
24  and  13  contained  in  configurations  (b)  have  40-degree  swept  wings  with  AR=2.1 
and  4  percent  thick  circular  airfoil  sections.  Configuration  13  is  similar  to  Configuration 
24  excepting  that  it  has  a  pressure  instrumented  store  instead  of  a  dummy  store  as  in 
Configuration  24.  Configuration  24  has  been  tested  at  the  free  stream  Mach  numbers 
of  0.925  and  0.950  at  0  degree,  2  degrees,  and  5  degrees  angles  of  attack.The  set  of 
experiments  on  Configuration  13  to  form  a  data  base  for  the  store  pressure  distribution 
computation  was  also  conducted  a',  the  .AEDC  4T  wind  tunnel  at  free  stream  transonic 
Mach  numbers  of  0.925,  0.952,  and  1.052  at  0  degree,  2  degrees  and  5  degrees  angles  of 
attack. 


In  the  earlier  analysis  we  obtained  good  agreement  between  theory  and  experiment 
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for  simple  wing-body  and  wing-body-pylon  configurations  as  could  be  observed  from 
Figures  9  to  14.  The  additional  complexities  like  wing  at  0  degTee  incidence,  thus  in¬ 
creasing  the  flow  accelerations  below  the  wing  which  is  the  region  of  interest  as  compared 
to  the  positive  angles  of  attack  of  the  wing  cases,  and  the  separated  store  in  the  close 
proximity  of  the  wing  and/or  pylon  involve  considerable  interference  effects.  Some  of 
these  are  first  order  and  some  others  are  second  order.  In  some  regions  of  the  flow-field 
first  order  theory  is  adequate,  whereas  second  order  corrections  are  required  for  some 
other  portions.  In  the  process  of  present  investigation,  we  identify  these  regimes  and 
indicate  the  ways  of  improving  the  results  further. 

AFFDL  configuration-24  with  wing- fusel  age- pylon  and  store  showed  considerable 
deviation  from  the  experiments  in  the  previous  analysis.  The  introduction  of  addi¬ 
tional  panels  along  with  their  better  distribution  improved  the  solution  for  the  flow-field 
velocities  upstream  of  the  wing  trailing  edge.  The  introduction  of  non-linear  effects  im¬ 
proved  the  results  slightly  (Figures  13  to  19).  Also,  in  this  case  the  equivalent  blowing 
method  to  account  for  weak  viscous-inviscid  interaction  outlined  in  Section  III  has  yielded 
the  results  as  presented  in  Figure  20.  The  strong  viscous  interaction  regions  near  the 
wing  trailing  edges  as  seen  in  Figure  20  are  beyond  the  scope  of  the  present  first  order 
theories  and  hence  the  sudden  accelerations  aft  of  x=18  and  the  subsequent  deceleration 
(Figures  14  to  19),  typical  of  the  near  wake  regions,  are  not  predicted  well.  The  thin  wing 
shapes  and  the  low  angles  of  attacks  involved  in  the  present  configurations  introduce  only 
minor  corrections  to  the  off  body  points  velocities  as  implied  in  equation  (29)  because  of 
the  rapid  decay  of  the  tangential  velocity  in  the  direction  normal  to  the  surface  of  the 
wing.  Also,  a  lower  Mach  number  of  0.82  was  used  in  the  linear  computation  instead 
of  the  actual  free-stream  value  of  0.925  to  examine  the  influence  of  local  compressibility 
coordinate  stretching  on  the  flow-field  velocities  for  linear  calculations.  These  results  are 
shown  in  Figure  21.  The  results  show  negligible  improvements  in  the  predicted  flow -field 
velocities.  Hence,  the  same  free  stream  Mach  number  as  the  wind  tunnel  test  value  has 
been  used  in  all  the  free  stream  subsonic  linear  computations.  The  transonic  non-linear 
corrections  introduced  by  the  integral  equation  method  as  described  in  Section  II  have 
been  incorporated  on  the  solution  obtained  after  the  viscous-inviscid  interaction  scheme. 
Thus,  we  can  observe  that  the  inclusion  of  weak  viscous-inviscid  interaction  together 
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with  nonlinear  corrections  improve  the  overall  prediction,  but  the  near  wake  regions 
could  not  be  adequately  simulated  with  the  present  scheme  as  expected. 

It  is  well  proven  that  the  panel  method  gives  good  results  for  subcritical  flow  over 
arbitrary  three-dimensional  geometries.  To  test  this  idea  again  for  the  present  case  of 
separated  store  configuration  a  trial  case,  configuration  001.  was  included  in  the  study 
at  a  free  stream  subsonic  Mach  number  of  0.7.  Figures  22  to  24  show  the  resubs  of  this 
calculation.  The  results  show  good  agreement  (Figures  22  and  23)  except  in  the  gap 
region  between  the  pylon  and  the  store.  The  two-dimensional  internal  flow  calculation 
as  described  in  Section  IV  has  been  subsequently  incorporated  for  the  flow  in  this  gap 
region  and  these  results  show  considerable  improvement  as  shown  in  Figure  24.  In  the 
absence  of  a  better  correlation  the  internal  and  external  effects  have  been  averaged  for 
the  present  cases. 

The  next  configuration  studied  incorporated  the  store  separation  effects  on  the  flow- 
field.  This  configuration  has  the  wing  at  2  degrees  angle  of  attack  and  hence  the  flow-field 
underneath  the  wing  does  not  experience  considerable  supercritical  accelerations.  Figures 
25  and  26  show  the  U  velocity  distribution  along  the  store  axis  for  various  positive  z 
locations.  Since  the  likely  internal  flow  in  the  gap  does  not  involve  appreciable  choking 
effects,  the  external  flow  computation  alone  yielded  good  predictions.  The  shock  location 
and  its  strength  were  not  determined  exactly.  Good  predictions  have  been  observed  at 
the  points  below  the  store  as  seen  in  Figure  27.  The  viscous  effects  were  considered 
insignificant  for  this  case  and  hence  have  not  been  implemented.  A  similar  study  to 
understand  the  influence  of  lowering  the  Mach  number  below  Mcr  for  linear  calculations 
yielded  slightly  different  solutions  as  shown  in  Figures  28. 

The  earlier  computation  results  on  configuration  221  which  showed  considerable 
deviation  between  theory  and  experiment  in  the  positive  i  values  were  identified  to  be 
due  to  the  presence  of  multiply  connected  regions.  This  error  was  fixed  by  introducing 
an  artificial  cut  connecting  the  wind  tunnel  wall  and  the  sting  store  junction.  The  results 
presented  in  Figures  29  and  30  show  significant  improvement.  The  artificial  cut  near  the 
pylon  trailing  edge  has  introduced  some  spurious  acceleration  aft  of  x=3.7.  The  linear 
and  non-linear  results  did  not  show  much  deviation  and  no  shock  was  observed  in  the 
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The  effects  of  the  store  are  more  significant  in  the  regions  between  the  store  axis  and 
the  wing  for  the  configurations  without  pylon  due  to  the  pressure  gradient  imposed  by 
the  wing  bottom  surface  in  addition  to  that  of  the  store  shape.  Hence,  the  configurations 
101  and  102  show  regions  of  shock  and  shock  boundary  layer  interaction.  In  addition, 
the  internal  flow  between  the  wing  and  the  store  becomes  an  important  factor  for  the 
flow-field  velocities  near  the  store  for  positive  i  values.  The  overall  prediction  for  Mach 
number  of  0.92  is  fairly  good  at  all  i  values  as  could  be  observed  from  Figures  31  to  36. 
Two  shocks, one  near  x=2.25  and  the  other  near  x=4.0,  are  closely  predicted.  Detailed 
comparisons  of  the  theoretical  and  experimental  V  and  W  velocities  are  avoided  in  most 
of  the  configurations  because  of  the  uncertainties  of  those  data  due  to  experimental 
intricacies.  When  the  Mach  number  was  increased  to  0.975,  the  flow  changed  its  behavior 
(Figures  37.38  and  47).  The  stronger  shock  that  occurs  near  x=4  makes  the  flow  fully- 
choked  and  hence  predominantly  internal.  The  shock  that  was  observed  near  x=2.25 
for  M=0.92  has  disappeared  in  this  case.  The  region  of  shock- induced  separation 
downstream  of  x==3.9  is  evident  from  the  experimental  data  given  in  Figure  47  from 
the  slow  recovery  of  the  pressure  downstream.  Although  the  present  calculations  show 
the  shock  near  x=4,  the  pressure  recovery  downstream  is  fairly  rapid  and  hence  the 
discrepancy  between  the  theory  and  experiment. 

The  influence  of  the  internal  flow  has  been  felt  in  the  inboard  and  outboard  sides 
of  the  store  as  well.  In  the  absence  of  any  correlation  between  the  correction  at  the 
store  centerline  and  at  other  regions  above  the  store,  no  correction  for  the  internal  flow 
computation  has  been  applied  at  these  locations.  Since  the  mass  flow  readjusts  itself  so 
that  the  throat  becomes  choked,  the  predicted  velocities  by  external  singularity  methods 
are  higher  than  the  date  as  shown  in  Figures  37  and  38.  The  source  correction  method 
used  by  other  investigators  to  simulate  the  flow  fields  near  the  engine  nacelle  could  be 
utilized  to  correct  the  inboard  and  outboard  stations.  Similar  technique  to  the  one  used 
for  configuration  101  has  been  applied  to  configuration  102.  A  thicker  store  in  this  case 
introduces  a  stronger  shock.  The  shock  location  was  upstream  of  the  actual  location  and 
its  strength  has  been  overpredicted  as  seen  in  Figure  39. 
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Figures  40  shows  the  results  of  the  introduction  of  pylon  to  the  previous 
configuration  along  a  section  outboard  of  the  axis.  The  non-linear  calculations  do  not 
show  the  shock  near  the  shoulder  of  the  store  which  might  be  due  to  the  fact  there  is 
no  sharp  gradient  of  the  linear  velocities  upstream.  Physically,  the  data  shows  a  separa¬ 
tion  bubble  on  the  stoTe  downstream  of  x=2.8  and  hence  there  is  an  acceleration  aft 
of  x=3.0.  Incipient  separation  downstream  of  x=3.25  is  observed.  These  details  could 
not  be  predicted  by  the  present  scheme  and  require  elaborate  computations  near  this 
section. 

Next,  the  results  of  configuration  121  at  Mach  numbers  0.92,  0.9T5,  and  1.025  are 
shown  in  Figures  41  to  46.  For  supersonic  free  stream,  consistent  with  what  has  been 
observed  in  the  shadowgraph  pictures  in  Reference  8  and  in  the  computed  results  from 
FL027  code,  a  detached  near  normal  shock  was  assumed  upstream  of  the  store  leading 
edge  and  the  adiabatic  shock  relation  there  yielded  a  Mach  number  downstream  of  the 
shock  to  be  0.9522.  The  non-linear  calculations  carried  out  at  this  Mach  number  did 
not  show  any  shock  or  significant  viscous  effects.  The  results  show  reasonable  agreement 
with  the  experimental  data. 

The  forward  and  reverse  sweep  technique  discussed  in  Section  rV  has  next  been 
applied  to  configurations  101  and  102  to  locate  the  shocks  properly.  Figures  47  and  48 
show  the  results  after  this  computation.  The  results  obtained  without  the  addition  of 
the  displacement  thickness  during  the  reverse  sweep  did  not  show  correct  shock  location 
for  Configuration  102.  The  addition  of  the  displacement  thickness  to  either  walls  indeed 
predicted  the  shock  location  fairly  accurately  for  both  the  stores.  Subsequent  study 
resolved  this  problem  by  directly  choosing  proper  artificial  viscosity  terms. 

The  results  for  the  store  pressure  computations  carried  out  using  the  method  dis¬ 
cussed  in  Section  V  are  next  presented  in  Figures  49  to  57.  It  is  interesting  to  note 
from  these  results  that  the  panel  method  alone  can  predict  the  store  on-body  pressures 
reasonably  well  in  regions  near  the  store  nose.  Also,  „nc  supersonic  pocket  near  the  nose 
shoulder  is  computed  well.  The  linear  theory  alone  fails  to  compare  well  with  the  experi¬ 
ments  due  to  the  interference  caused  by  the  presence  of  the  pylon  and  wing.  Whereas  the 
present  theory  does  a  good  job  of  predicting  the  two  accelerative  and  decelerative  behaviors 
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in  the  center  and  aft  portions  of  the  store.  Overall  favorable  agreement  between  the 
present  computational  scheme  and  the  experimental  results  has  been  observed  along  the 
bottom  periphery  of  the  store.  As  the  free  stream  Mach  number  has  been  increased  from 
free  stream  subsonic  to  supersonic  values,  the  acceleration  flow  near  the  nose  shoulder 
apparently  shows  a  strong  shock  which  has  not  been  predicted  well  in  the  linear  com¬ 
putations  and  hence  does  not  show  up  in  the  present  computations  as  well.  The  matching 
of  the  linear  velocities  alone  up  to  the  store  nose  shoulder,  i.e.  up  tox=2.0  and  the  present 
solution  of  adding  the  interference  flow  to  the  full  potential  solution  could  be  a  viable 
approach  to  obtaion  a  better  approximation  for  the  present  problem.  As  mentioned  in 
Section  V,  the  linear  addition  of  the  influences  due  to  the  non-linear  effects  (both  due 
to  the  governing  differential  equation  and  due  to  the  interference  effects)  is  not  fully 
justifiable.  However,  an  attempt  has  been  herein  made  to  try  such  naive  methods  to  un¬ 
derstand  the  physical  significance  of  each  of  the  nonlinearities.  The  results  presented  in 
Figures  49  to  57  indicate  that  the  above  method  proves  that  the  present  method  is  good 
for  a  quick  and  a  fairly  accurate  estimation  of  the  store  pressure  distribution  and  a  more 
refined  approach  which  builds  up  a  correlation  between  the  linear  and  non-linear  effects 
independently  for  the  two  parts  and  then  implementing  a  cross  correlation  between  the 
two  should  yield  more  accurate  results.  This  methodology  is  under  present  investigation. 
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SECTION  vn 


CONCLUSIONS  AND  RECOMMENDATIONS 


Conclusions 

From  the  above  analysis,  the  following  conclusions  are  appropriate: 

(1)  The  integral  equation  method  of  transonic  non-linearity  correction  applied  over  the 
panel  method  linear  velocities  is  adequate  for  the  flow-field  predictions  at  points 
free  of  very  pronounced  transonic  interference  effects. 

(2)  The  weak  viscous  effects  on  flow  field  velocities  for  thin  wings  as  in  the  present 
configurations  are  small  and  hence  can  be  ignored. 

(3)  Significant  non-linear  effects  caused  both  due  to  complex  geometry  and  due  to 
transonic  flow  non-linearities  can  be  locally  corrected,  as  discussed  in  this  report, 
to  yield  approximate  results.  This  approach  avoids  extensive  computing  that  is 
otherwise  necessary  by  a  more  rigorous  procedure. 

(4)  Efforts  to  correct  the  store  on-body  pressure  distribution  by  a  simple  procedure, 
as  discussed  in  Section  V,  did  not  yield  very  good  predictions.  A  more  elaborate 
empirical  approach  is  deemed  necessary. 

Recommendations 

(1)  The  incorporation  of  a  procedure  for  exact  wake  shape  determination  is  needed. 
This  could  be  done  by  iterating  approximate  wake  shapes  successively  using  the 
boundary  conditions  of  tangential  velocity  continuity  and  normal  velocity  discon¬ 
tinuity  across  the  wake  by  using  the  same  singularity  strengths  on  the  wing  for  each 
iteration. 

(2)  Implementation  of  the  weak  separation  effects  into  the  panel  method  by  the  ap¬ 
propriate  modification  of  the  present  viscous-inviscid  interaction  code  is  needed. 

(3)  MacCormack’s  recent  implicit-explicit  finite  difference  scheme  needs  to  be  imple¬ 
mented  into  the  VNAP  code  of  internal  flow  computation  to  reduce  total  computing 
time. 

(4)  Development  of  a  good  correlation  between  the  center  section  velocity  and  velocities 
at  other  points  around  the  store  is  needed  in  order  to  correct  the  velocities  at  the 
latter  points  in  the  internal  flow  computation. 


(5)  Finally,  the  various  individual  codes  discussed  in  the  different  sections  need  to  be 
combined  to  make  a  unified  computational  scheme  with  necessary  user  options. 
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Figure  3.  Modified  Airfoil  Geometry  with  Separation  Effects  Correction 
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Figure  9.  U  VELOCITIES  FOR  M  =  0.925  ;  ALPHA  =  0  dog 
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Figure  13.  V  VELOCITIES  FOR  M  »  0.925  ;  ALPHA  =  0  deg. 


48 


CONFIGURATION-24 
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Figure  20.  U  Velocities  at  Y-4.2S,  Z=- 1.23  for  Configuration  !■ 

M=0 . 925 ,  a=5° 
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Figure  25.  U  Velocities  at  Y=0.0,  Z=0.612  for  Configuration  201 

M=0 . 92 ,  j:=2° 
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Figure  42.  U  Velocities  in  the  Gap  at  Y=0.0,  Z=0.392  for  Configuration  121 

>1=0.975,  7=0° 


77 


A  A  A 


>1=1.025 


□  cl 

©  UT 

A  UE 


0.0 


§ 

A 


g 


8 


a 

0 


0 


A 

§ 


A 


1.0 


2.0 


3.0 


5.0 


Figure  46.  U  Velocities  at  Y=-0.35,  Z=0.392  for  Configuration  121 

M= 1.025,  ;=0° 
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Figure  47.  V  Velocities  in  the  Gap  at  Y-0.0,  Z=0.39:  for  Configuration  101 

>1=0.97  5,  J=0° 

(With  or  Without  Displacement  Effect) 
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Surface  Along  the  Axis,  M=0.S52,  a-0 


Figure  52 •  C  on  th6  Store 
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Figure  57.  C  on  the  Store  Surface  Along  the  Axis,  M= 1.052,  3=5 
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All  dimensions  in  inches 
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Detail  and  Dimensions  of  the  Stores  Used 
Configurations  (a)  (continued) 
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WIND  TUNNEL  SIDEWALL 
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ALL  DIMENSIONS  IN  INCHES 
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Wing  Details  and  Dimensions 
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Pvxon  Details  and  Dimensions 
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6ose  Coordinates 

4-percent 

Airfoil 

X ,  in. 

R ,  in. 
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t/c.  oercent 
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APPENDIX  B 

FORTRAN  LISTING  OF  THE  NON-LINEAR  CORRECTION  PROGRAM 
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THIS  IS  THE  PROGARM  TO  INCLUDE  THE  NON-LINEAR  EFFECTS  TO  THE 
LINEAR  VELOCITIES  OBTAINED  FROM  THE  PANEL  METHOD. 


IMPLICIT  REAL*8!A-H.O-Z> 

DIMENSION  XDUMMY I  2013)  .  YDUMMVI  300)  , ZDUMMY(300>  ,CPDUM(300) 

DIMENSION  UOUMMY I  300 ) .VOUMMY!  300) .WDUMMY! 300) 

DIMENSION  X<600> . Y ( 600 ) . Z ! S00 > , U ( 600 > , VI 600) ,W( 600) .CPI 600) 

DIMENSION  WEAR  I  600  > 

DIMENSION  UBARI 600) , YBAR ! 600) .ZBAR ( 600) 

DIMENSION  VBARI600) 

DIMENSION  VSQL 1 1  600  > , VSQNL I ( 600 ) ,CPLI ( 600 ) , CPNL I ( 600 ) 

COMMON  /POINT/  A.B.C 
COMMON  /FIELD/  UDUMMY . VOUMMY .WDUMMY 
COMMON  /VEL/UL.II 

******  M  IS  THE  NUMBER  OF  FIELD  POINTS  FOR  WHICH  THE  NONLINEAR  CORRECTION 

******  IS  SOUGHT. 

******  N  IS  THE  NUMBER  OF  DUMMY  POINTS  USED  IN  I NTEGRATION . PRESENT 

******  CASE  IT  IS  6*6*6  *  216  POINTS 
******  aMACH  IS  THE  FREE  STREAM  MACH  NUMBER 
******  NSTO  IS  FOR  THE  STORE  INFORMATION. 

******  NSTO  EO  1  WHEN  STORE  IS  PRESENT 

******  IF  NST0  Is  1 . THEN  XSTO.YSTO.ZSTO  HAS  TO  BE  GIVEN  IN  ORDER 
******  T0  CONVERT  THE  COORDINATES  IN  THE  OUTPUT  FROM  WING 
******  COORDINATE  SYSTEM  TO  STORE  COORDINATE  SYSTEM 

******  IF  NSTO  IS  NOT  I , THEN  XSTO . YSTO. ZSTO  DOES  NOT  HAVE  TO  BE  GIVEN. 

******  VELOC  IS  THE  FREE  STREAM  VELOCITY  IN  FT/SEC 

******  EITHER  THIS  CAN  BE  UNITY  IN.  WHICH  CASE  THE  OUTPUT  VELOCITIES  ARE 
......  U/U  INFINITY.  OR  IT  CAN  BE  THE  ACTUAL  VELOCITY  AND  IN  THIS  CASE 

......  THE  OUTPUT  VELOCITY  WILL  BE  IN  FEET/SEC  OR  IN  THE  GIVEN 

......  UNIT  AS  DEFINED  BY  VELOC 

******  U.V.W  ARE  THE  PERTURBATION  VELOCITY  COMPONENTS  IN  PHYSICAL  COORDINATES 
......  UBAR  , VBAR , WBAR  ARE  THE  REDUCED  VELOCITY  COMPONENTS 

******  X.Y.Z  ARE  THE  POINTS  IN  PHYSICAL  COORD IATE  FOR  WHICH  THE  VELOCITY 
......  DISTRIBUTION  SOUGHT. 

******  XBAR, YBAR. ZBAR  ARE  THE  REDUCED  COOR D I  NATES . FOR  REFERENCE  SEE  NACA 
******  REPORT  1217 

..*...  UDUMMY  .VOUMMY  .WOUMMY  ARE  THE  PERTURBATION  VELOCITY  COMPONENTS 
******  AT  THE  DUMMY  POINTS  IN  THE  SCALED  SPACE  FOR  .INFINITY  TO  -INFINITY 
......  ATTACH  THE  OUTPUT  OF  HESSPROGRAM  FOR  THE  DUMMY  POINT  AND  FOR  THE 

......  field  POINTS  AFTER  THE  FIRST  OATA  IN  THIS  PROGRAM 

CALL  ASSIGN  <5, ’DB0:NEW240NL . IN' ) 

CALL  ASSIGN  17. ’  DB0: NEW240NL . OUT ' ) 

READ  ( 5. 345  )  M.  AMACH, VELOC, NSTO. ALPSUN 
OATA  M, AMACH. VELOC . NSTO . AL PSUN/ 86 , 0 . 925 , l . , 1 ,0.0/ 

DATA  M, AMACH, VELOC, NSTO/ 86. 0.92 5. 1.0,0./ 

IF  (  NSTO.EQ.l)  READ  ( 5 , 55 ) XSTO , YSTO , ZSTO 
XSTO-5 . 

YSTO-6. 

ZSTO*-  1  . 

55  FORMAT  (3F10.6) 

N  -  226 

•  345  FORMAT! 1 15 , 2F 10.6. 15, F 10.6) 

READ  (5,81 HXII ) ,Y< I > ,Z( I ) ,U( I ) ,V< I > ,W( I >,CP( I) .1-1 ,M) 

81  FORMAT! 7F 10. 6) 
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READ  <5,81 ) <  XDUMMY ( I ) , YDUMMY C I > .ZDUMMY! I ) .UDUMMY! I ) . 

1  VDUMMY  (  I  )  , UDUMMY <  1  )  .CPOUMC  1)  ,  1-1  ,  N> 

C  WRITE  <6,342) 

342  FORMAT!  5X,  ’  X  ’  ,  1 IX ,  '  Y '  ,  1 IX ,  ’  Z  ’ ,  10X ,  ’  UL  ’  ,  6X  , 

1  6X,’U’,12X,’VL’,11X,'V’,11X,’WL',11X,’W,,7X, ’CPL1N  * ,3X  ,  ‘CPN  '  ) 

CONST “ (AMACH*AMACH  >*2.4/1  .0 
C  CONST“AMACH“l . 1*2. 4/1  .0 

ANX-DCOS!ALPSUN*3. 14159265/180.  i 
ANZ“DSIN!ALPSUN*3. 14159265/180. > 

WRITE  (7,3285) 

3285  F  ORMAT  <  5X , ' THE  CONFIGURATION  IS  24  .0.925’) 

DO  3967  1*1,216 

UDUMMY ( I ) “UDUMMY  < I )-ANX 

VDUMMY! I ) “VDUMMY ( I ) 

WDUMMY ( I ) -WDUMMY  < I ) -AN2 
3967  CONTINUE 

DO  840  I-l.M 
U  < I ) “U ( I ) -ANX 
V  < I  )  «V  < I > 

W< I >-W( I >-ANZ 
840  CONTINUE 

BETA=DSQRT ( 1 . 0-AMACH*AMACH > 

DO  2034  I -  1 ,216 

C  UDUMMY ( I ) “UDUMMY ( 1 )*CONST/(B£TA*BETA> 

UDUMMY! I ) “UDUMMY ( I )*CONST 
2034  CONTINUE 

DO  10  I “ 1 , M 
YBAR! I >»Y< I ) ‘BETA 
ZBAR! I  )  -Z ( I  ) ‘BETA 

C  UBAR ( I )-U( I )*CONST /(BETA* BETA) 

UBAR ( I > -U ( I ) ‘CONST 
C  PRINT  3490, U< I ) .UBAR! 1 > 

3490  F  ORMAT ( 2F 1 5 . 6 ) 

C  VBAR ( I ) - V! I ) ‘CONST /(8£TA**3) 

VBAR! I > “V ( 1 ) ‘CONST /BETA 
C  WBAR ! I )-W{ I )*CONST/( 8 ETA** 3 ) 

WBAR ( I l-WBAR! I ) ‘CONST/BETA 
A-X! 1 ) 

B-YBAR! I ) 

C-ZBAR! I > 

D-UBAR!  I  ) 

CALL  TRIPLE(E) 

11*1 

G-UBAR! I )♦( 0*0) /2.D0-! 0.079577 )‘E 
G 1  - <  2 . 0*U8AR ( r )-l .0) 

G2- < (2.0*0.079577‘E)-(G1 )) 

IF (UBAR! I ) .GT. 1 .0)  G“( 1 .0*DSQRT< ( 2 .0*0 ,079577‘E ) - ( G 1 1 > ) 

IF ( UBAR ( I ) .GT. 1 .0.AND.G2.GT.0.0)  G-! 1 .0+DSQRT!G2> > 

IF ( UBAR! X  > . L T . 1 . 0 . ANO . G2 . GT . 0 . 0 )  G“( 1 .0-DSQRT(G2  > ) 

IF (UBAR! D.LT.1.0)  G“(l . 0-DSQRT ( ( 2 . 0«0 . 079577*E ) - ( G 1 ) > ) 

IF (UBAR! I ) .EQ. 1 .0)  G-UBAR(I) 

IF (G2 . EQ.0.0)  G-1.0 

I F ( G2 . L  T . 0. 0 )  G-UBAR! I )+!D*D>/2. 0D0-! 0.079577 >*E 
WA-DSQRT ( 3 . D0 ) 

G* (WA  + 1 .0)*! 1 .0-! 1 . 0- < WA- 1.0)*UBAR!  I ) >  “0 . 5 > 

FORMAT!5X.2F15.8) 

UL-G 

CALL  TRIPLE  !E> 

II-II+l 

IF! 1 1 .GE.3)  GO  TO  1001 
GO  TO  1000 
1001  CONTINUE 

CALL  TRIP(SUM) 

VFIN-VBAR! I ) -0.079577*SUM 
CALL  TR I  PE ! TOT ) 

WF IN-WBAR! I . -0 . 079677«TOT 
C  TRANU“(G*BETA*BETA> /CONST 

TRANU-G /CONST 

C  TRANV«(VFIN*BETA“3)/CONST 


C 

1000 

C 


c 

c 


c 

c 

1005 
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TRANV-VFIN*BETA/CONST 
C  TRANW- (Wf 1 N*BETA**3 ) /CONST 

TRANW-WF IN*BETA/CONST 

TRANU*  < TRANU  +ANX  > -VEL  OC 
TRANV* ( 0 . 0+ TRANV  >  *VE  LOC 
TRANW- ( AN Z+ TRANW ) * VE  L  OC 
U<  I  >  *  ! ANX  +  U  <  I  >  )*VElOC 
V( I  )-V(  I ) * VELOC 
W ( I } * ( ANZ+W! I ) ) -VELOC 
V( I )»V( I > 

TRANV-TRANV 

VSQL I ( !)-U< I )**2+V( I >**2+W( I >**2 
VSQNL I ( I )-TRANU**2+TRANV**2*TRANW**2 
IF  (NSTO.NE. 1 )  GO  TO  100 
X ( I )  *X ( I  )  -XSTO 
Y< I ) -- YSTO+Y  < I ) 

2  < I ) -Z  < I ) -ZSTO 

C100  WRITE  <6,234)X< I ) , Y< I ) ,Z( I > ,U< I > . TRANU, V( I > , TRANV, W 

100  CONTINUE 

C  1  TRANW 

C  WRITE  (6.3450JE 

3450  FORMAT! 10X.F15. 6) 

WRITE  ( 7.231 >X! I ) . Y! I ) ,2! I >,U< I >, TRANU. V! 1 > .TRANV.W! I 
231  FORMAT ( 6F 10. 6 } 

234  F  OR MAT !4F12.6,5F12.6) 

10  CONTINUE 

STOP 
END 

SUBROUTINE  TR I  RLE < RESUL T > 

C 

IMPLICIT  REAL  *8  <A-H.O-Z>  ^ 

DIMENSION  P<8> ,H!8) 

DIMENSION  UDUMMY ( 300) , VDUMMYl 300 >, WDUMMY < 300 > 

COMMON  /POINT/  X0.Y0.Z0 
COMMON  /FIELD/  UOUMMY . VDUMMY . WDUMMY 
COMMON  /VEL/  UL , I  I 
C 

P< 1 1*0.93246951420315200 
P ( 2  )  *-P (  1  ) 

P ! 3 ) *0 . 66 1 20938646626500 

p { 4  )  »-P ( 3  > 

P(  5 1-0. 23861 91 86083 197D0 
P'6)--P(5> 

H ( 1 ) *0 , 17132449237917000 
H ( 2  >  *H I  1  ) 

H( 3 1*0. 3607615730481 3900 
H ( 4 ) *H ( 3  > 

H! 5 ) *0. 46791 39345 72691 D0 
H ( 6  >  *H ( 5 ) 

RESULT-0. D0 

DO  500  1*1,6 

DO  500  0-1,6 

DO  500  M-l ,6 

K-< ( 1-1 >*36+10-1 >  *6  +M ) 

IF! II .EQ.  1 ) UL -UOUMMY ! X  > 

RESULT-RESULT+H! I > -H ( 0  > *H ( M > *F < P ( I > . PI  0 > , P ( M ) > 

500  CONTINUE 
RETURN 
END 

DOUBLE  PRECISION  FUNCTION  FiX.Y.Z) 

C 

IMPLICIT  REAL  *8  (A-H.O-Z) 

COMMON  /POINT/  X0.Y0.Z0 
COMMON  /VEL/  UL  .  I  I 
C 

PI-3. 14159265400 

PII-PI/2 .000 

A* ( X0-OTAN ( P1I*X)>**2 


ID. 


>. TRANW 


■j  v-Srim..-  - 
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B  -  ( Y0-DTAN ( P 1 1  * V  > )  **2 
C-IZ0-OTANIPI I*Z>  >**2 
D-I 1 .0/DCOSIPI I*X> )  **2 
E-I 1 .0/DCOSIPI I*Y> >  *  *2 
G-I 1 .0/OCOSI P I  I  *Z  > 1**2 

F-((PI**3>/8.  >*!<UL**2>/2.>*( ( 2 .0*A ) -B-C > * ( D*E*G ) 

1  / I  I  A*  B+C ) *  *2 . 5 ) 

RETURN 

END 

SUBROUTINE  TRIPIRESULT) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  P ( 8  > , H <  8  > 

DIMENSION  UBAR1600) 

COMMON  /POINT/  X0.V0.Z0 
COMMON/VEL/UL . I  I 

DIMENSION  U DUMMY ( 3001 . VOUMHYC  300 J .VOUMM Y(300) 

COMMON  /FIELD/  UOUMMY , VDUMMY .WDUMMY 
PI  1 >*0.93246  9514  20315200 
P(2)*-P( 1 > 

P (3 1=0.66120938646626500 
P  I  4  )  *  -P  I  3  > 

PI  5  )-0. 2366191860831 9700 
P I  6  >*-P I  5  > 

HI  1  >*0.  1  71 324  4923791 70D0 
HI21-HI 1 > 

HI  3 >*0.3607  6157304813900 
HI  4  >-Ht  3  > 

HI  5 >-0.46791 393457269100 
H I  6 ) *H I  5 ) 

RESULT-0.00 

DO  500  1-1.6 

DO  500  J-l .6 

DO  500  M-l .6 

K-(t I-l>*36+< J-l >*6+M> 

UL  -UOUMMY  I  K  > 

RESULT-RESULT-HI I>«HIJ>*H(M>*FF(P(I).P(J>.P(M>> 

500  CONTINUE 
RETURN 
END 

DOUBLE  PRECISION  FUNCTION  FF1X.Y.Z) 

IMPLICIT  REAL *8  (A-H.O-Z) 

COMMON  /POINT/  X0.Y0.Z0 
COMMON  /VEL/  UL  ,  1 1 

PI-3. 141592654D0 
PII-PI/2.000 
A- I X0-OTAN I P I I*X> 1**2 
B » 1 Y0-DTAN I P I  I  * Y  > >**2 
C*(Z0-OTAN( P 1 1 -Z ) >**2 
D-I 1 .0/DCOSI P I I*X  >  >**2 
E  * ( 1 .0/DCOSIP 1 1 *Y  > 1**2 
G-I 1  ,0/DCOS(PII*Z) )**2 

FF-I (PI**3>/8. )  *  I IUL**2>/2. >* ( 3 . 0*DSQRT ( A ) -OSOR T( 8 > >*<0*E*G> 
/( I A-B+C  >  **2 . 5 ) 


RETURN 
END 

SUBROUTINE 


TR  I  PE(  RESULT  > 


IMPLICIT  REAL *8  IA-H.O-Z) 

DIMENSION  P I  8 ) , H ( 8  > 

DIMENSION  U8AR ( 600 1 

COMMON  /POINT/  X0.Y0.Z0 
COMMON  /VEL/  UL  ,  II 

DIMENSION  U DUMMY ( 300 ) .VDUMMY ( 300 ) . WDUMMY ( 300) 
COMMON  /FIELD/  UOUMMY , VDUMMY , WDUMMY 


PI  1  ) -0.93246951 4203 1 52D0 
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_  -i 


P(2)»-P<  1  > 

PI  31*0.66120938646626500 
P ( 4 ) »-P ( 3 ) 

P( 5) -0.23861 9 1860831 9700 
P(6>*-P15> 

H! 1 1-0. 1713244923791 70D0 
H ( 2 ) *H  < 1  ) 

H<  3  )-0. 360761 573048 139D0 
H ( 4 ) *H ( 3 1 

H(S  >-0. 4  6 791 39345 7269 1D0 

H ( 6  >  *H( 5 ) 

RESULT-0.00 

00  500  1*1,6 

DO  500  J-l . 6 

DO  500  M*1  ,6 

K-< ( I-l >*36*<U-1 >*6«M> 

C  UL-UDUMMY(K) 

RESULT-RESULTfH! I>»H<0>*H<M)*FFF<PU).P<J>,P<M>) 

500  CONTINUE 
RETURN 
END 

DOUBLE  PRECISION  FUNCTION  FFFIX.Y.Z) 

C 

IMPLICIT  REAL  *8  <A-H,0-Z> 

COMMON  /POINT/  X0.Y0.Z0 
COMMON  /VEL/  UL , I  1 
C 

PI-3. 1 4 1 592654D0 
PII-PI/2.0D0 
A* ( X0-DTAN ' P 1 I*X  > 1**2 
B*(Y0-DTAN(PII*Y) 1**2 
C*(Z0-DTAN(PII*Z) )**2 
D* ( 1 . 0/DCOS ! P I I*X) 1**2 
E * ( 1 . 0/DCOS I P I  I *Y ) ) **2 
G-< 1 .0/DCOS! PI  1*2) 1**2 

FFF*((PI**3)/8. )«! <Ul**2>/2. 1 • < 3 .0*OSORT< A > «DSORT< C 1 >«<D*E*C> 
1  /( (A»B*C)**2.5) 

RETURN 

END 
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THE  FOLLOWING  TWO  SUBROUTINES  CALCULATES  THE  PRESSURE  COEFFICIENTS 
ON  AND  OFF  THE  BODY  BASED  ON  THE  SLENDER  BODY  ASSUMPTIONS. 


SUBROUTINE  CPOF F ! XC , YC , ZC . VXA . VYA , VZA > 

DIMENSION  ALP HAX!  10) .ALPHAY! 10) , AL  P HAZ ( 10) 

COMMON  /OUTPUT/FPR INT, IANGLE 

COMMON  /ANGLE/  MANGL E , AL P HAX . AL PHAY . AL P HAZ , SUNBET 
COMMON  /  SCALE  /  FC.  BL 
AMACH-SQRT! 1 . -SUNBET'SUN BET ! 

SUNBET  -SQR  T (  1 . -AMACH*AMACH ) 

UA-VXA-ALPHAXI IANGLE ) 

VA  =.V  YA 

WA- VZA-AL  P  HAZ ( IANGLE ) 

VX -COS ! ACOS  <  AL  P  HAX  < IANGLE ) > /SUNBET >♦<  UA/ <  SUNBET'SUNBET > ) 
VY-VA/SUN8ET 

VZ-SI N<ASIN<ALPHAZ< IANGLE ) ) /SUNBET >*< WA/SUNBET) 

AG 1-UA/! SUNBET'SUNBET ) 

AG2-VA/SUNBET 
AG3-WA/SUNBET 
VSQ-VX"2*VV"2*VZ"2 
IF (SUN8ET. EQ. 1 . )  GO  TO  3290 

CP  =*  (  1  .  /(0. 7*AMACH*AMACH  )  )*  (  (  1  .  + ! AMACH'AMACH/ 5  .  >*!  1  .  -VSQ)  >"3. 5 
1  -1  .  ) 

GO  TO  3300 
3290  CP-1. -VSQ 

3300  YC-YC/SUNBET*FC 

ZC-ZC/SUNBET*FC 
XC-XC'FC 

WR ITE  <  7 . 540)  XC .YC .ZC . VX , VY . VZ .CP 
540  FORMAT! 7( 2X , F 10.6 ) ) 

RETURN 

END 

SUBROUTINE  CPSUNO ( XC . YC . ZC . VXA . VYA . VZA ) 

DIMENSION  ALP HAX!  10) , ALP  HAY!  10) .ALPHAZ! 10) 

COMMON  /ANGLE/  MANGLE .ALPHAX .ALPHAY. ALPHAZ .SUNBET 
COMMON  /OUTPUT/  F PR  I  NT . I ANGL E 
COMMON  /  SCALE  /  FC.  BL 
AMACH*SQRT ( 1 . -SUNBET'SUNBET ! 

UA-VXA-ALPHAX! IANGLE ) 

VA-VYA 

WA-VZA-ALPHAZ! IANGLE) 

V:<-COS(ACOS!ALPHAX<  IANGLE)  ) /SUNBET  )*(  UA/ f  SUNBET'SUNBET  )  ) 
VY-VA/SUNBET 

VZ -SIN! ASIN! ALPHAZ! IANGLE) ) /SUNBET )*< WA/ SUNBET > 
AGI-UA/(SUNBET*SUNBET) 

AG2-VA/SUNBET 
AG3- WA/SUNBET 
VSQ-VX**2*VY**2*VZ**2 
IF (SUNBET. EO. 1 . )  GO  TO  3290 

CP  * !  1  . / ( 0 . 7*AMACH*AMACH ) )*! ( 1 . ♦ !AMACH*AMACH/5 .  )*! 1 . -VSQ) >**3.5-1 
1  .  ) 

GO  TO  3300 
3290  CP-1. -VSQ 

3300  YC-YC/SUNBET*FC 

ZC-ZC/SUNBET*FC 
XC-XC  *F  C 

WRITE! 7.345)  XC . YC . ZC . VX . VY . VZ , CP 
345  FORMAT! 5X . 7! F 10. 6 ,2X ) ) 

RETURN 

ENO 
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THIS  PROGRAM  READS  AND  REARRANGES  THE  INPUTS  FROM  THE  CIRCUMFERENTIAL  TO 
AXIAL  ORDER  FOR  NON-LIFTING  SECTIONS.  NOTE  THAT  A  MAXIMUM  OF  THREE 

NON-LIFTING  BODIES 

COULD  BE  PROCESSED.  NOTE  ALSO  THAT  A  CORRECTION  FOR  THE  WAKE 
THICKNESS  EFFECT  HAS  BEEN  INCORPORATED  BY  ADDING  6*H  AS 
WAKE  LENGTH,  WHERE  H  IS  THE  SEMI  BASE  HEIGHT. 

LOGICAL  LASSTR I . NEWSTR 1 .OFF , SECT , LASS EC . F USE . WAL L .EVEN 
DIMENSION  XL  ASTI  (  100)  ,  YLAST1  (  1 J0TJ0T  >  .ZLAST1  (  100)  , 

1  X  LAST  IP  1 ( 100 ) . YLAST1 P 1 ( 100!  , Z LAST  IP  l  (  130) , 

2  XLAST I P  2 (  100) . YLASTIP2* 100) ,ZLAST1P2<  I  00 ) , 

3  XOLD*  20 . 80  > .VOID! 20. B0).2OLD(20,80).MLIN1(100), 

4  ML  I N 1 P ( 100 ) , XNEW1 B (  10, 100 ) , YNEW1 B  <  10. 100) , 

5  ZNEW1B* 10, 100! .XNEW1T* 10, 100) .YNEW1T* 10. 100) . 

6  ZNEW1T* 10, 100) .LNEW1B* 10. 100) .MNEW1B* 10. 100) . 

7  LNEW1T* 10 , 100) . MNEW1T* 10 , 100 ) .XNEWl ( 1250) . XNEW2* 1250) , 

8  YNEW1 < 1250) , YNEW2 (  I  2 50 )  . ZNEW1 (  1 250 ) , ZNEW2*  1250), 

9  L NEW  11  1250) . L  NEW2 (  1 250 )  .MNEW1 (  1  250 )  , MNEW2 (  1250)  . 

A  X1I650) ,Y1 *650) ,Z1 *650) ,X2( 650) ,V2<650) ,Z2<  650) . 

B  L 1 ( 650 > ,  Ml(650),L2(650),MZ(650> .TITLE* 17) , 

C  X(  1250)  ,Y( 1250)  ,Z(  1250) . L  <  1250)  ,M(  12501 . 

D  NSORCE ( 10) , NWAKE (  10)  , NSTR IP (  10) . IPCV!  10)  . I  X  F  LAG '  10) .AL< 10)  . 

E  IG1 ( 50, 10) , IGN( 50,  10 ) .WIDTH! 50, 10) . WIOKTP ( 2 . I  0 >  . NXT( 50 > , 

F  NSPL IT< 10) ,NL I «E1  ( 10) . NL INENI 10) 

ASSIGN  THE  UNITS  TO  READ  THE  DATA. 

CALL  ASSIGN* 1 . •DBA0:TEST1 . IN' > 

CALL  ASS  I GN ( 3 , ’ OBA0 : TEST  I . OUT '  ) 

CALL  ASSIGN<5, 'DBA0iTEMP.DAT' ) 

TYPE*. 'PLEASE  SAY  I  TO  INOICATE  THAT  BL.FC.RE  VALUES  ARE  INPUT' 

ACCEPT* , ISTOP 

IF* ISTOP .NE. I )  STOP 

FIRST  READ  AND  WRITE  OOWN  THE  TITLE  AND  THE  CONTROL  FLAGS. 


READ* 1.1)  TITLE 
WR I TE ( 3 , 1 )  TITLE 
FORMAT* 1 7A4 I 

READ* 1,1234)  CASE . L I F S EC , MOMENT , KUTTA , NOFF  , L 1ST ,MPR . I  OUT, IG .LASWAK, 
I A TACK , IWIDTH, ISAVE.SYM1 , S YM2 , AMACH . F C , BL 
WR I T  E  < 3 .2  )  CASE , L IF  SEC .MOMENT .KUTTA, NOFF .L 1ST  ,MPR ,  10 UT. IG .LASWAK, 
IATACK ,  IWIDTH. I  SAVE .SYM1 ,SYM2 . BL . FC .AMACH 
FORMAT* A4 ,I2I3,3X,3F4.2,F12.3,F5.3) 

F  ORMAT  *  A4 , 12I3.3F4.0.F12.0.F4.0) 

READ<1.3)  IAL <  I  )  . I  -  1  , IATACK ) 

WR I TE ( 3 , 3 )  <AL ( I ) . I *1 , IATACK) 

FORMAT*  6F 10.6) 

I F ( MOMENT . NE . 0 )  THEN 


no 


READ (1.4)  0RIGNX,0R1GNY,0RIGNZ 
WRITE13.4)  ORIGNX.ORIGNY.ORIGNZ 
ELSE 
END  IF 

4  FORMATI3F10.6) 

READ!  1  , 1235  >  < NSORCE! I i  , NWAKE! I ) . NSTR I P < I > . NSPL I T( I) . NL I NE l ( I > , 

1  NL I NEN (  I  ) , I XF  LAG ( I > . I PCV(  I > , I ■ 1  .LIFSEC) 

UR  I TE  < 3 . 5 )  ( NSORCEI  I > , NWAKE I  I ) . NSTR IP (  I > , IPCVI I )  , I XF  L AG  I  I ) , I - 1 , L I F SEC ) 

5  FORMAT! 41 31 4 , 21 2  )  ! 

1235  FORMAT! 21 7  I  4 . 12 ) ) 

IF! I G . L  E .  <n  GO  TO  6 
DO  7  0-1 .LIFSEC 
0STRIP-NSTR1P! J ) 

READ!  1,8)  <IG1!I,J>.IGN!I.0>.I-1  . JSTRIP  ) 

7  WRITEI3.8)  !  IG1 (  I , J ) . IGNI I ,0 ) . I  =  1 . JSTRIP  ) 

8  F  OR MAT ! 1614) 

6  IF C IWIDTH . EO.0)  GO  TO  9999 
DO  800  J= 1 , LIFSEC 

JKK-0 

XX  -  I XF  LAG ! J ) 

IFIXX.EQ.0)  GO  TO  9 
IF!KX-2>  1200,1300,1200 
1200  OKX-1 
GO  TO  9 

1300  JKX-2 

9  OS-OSTR I P -0  KX 
OP  1 -JS  + 1 

READ! 1,3)  UIDKTR! 1.0) , (WIDTH! < , J ) , K-2 , JP1 ) .WIDKTR! 2.0 ) 

WRITE (3, 3)  WIOKTR! 1 , 0 ) , ( WI DTH ( K . 0  > . K -2 . OP  1  > , WI DKTR (2,0) 

800  CONTINUE 

9999  CONTINUE 

C 

C  NEXT  READ  THE  PANEL  AND  OFF  BODY  COORDINATES. 

C 

1-0 

7654  1-1*1 

READ! 1 , 1001X1! I ) ,Y1( I) ,Z1( I ) ,L1(I ) .Ml ( I ) ,XZ< I ) , Y2< I) . 

1  Z2( 1 ) ,L2( I ) ,M2< 1 ) 

IF (L  1  (  I  ) . EO. 3)  THEN 
NONBDY-! 1-1 )*2*1 
ION-  1 

GO  TO  8765 
ELSE 

IF ( L2< I ) . EQ. 3 )  THEN 

NONBDY-I-2 

ION-  I 

GO  TO  8765 
ELSE 

GO  TO  7654 
END  IF 
ENO  IF 

8765  0-1 

5432  0-0*1 

READ! 1,201)  X1(0),Y1(0),Z1(J>,H(J).X2(0).Y2(J),Z2(0),L2(J) 

I F ( L 1 (0 ) . EO . 3 )  THEN 
NOF  BOY- ( 0  -  I -  1 )*2*1 
I  OF  -0 

GO  TO  6543 
ELSE 

IF(L2(J)  .EQ.3)  THEN 
NOF  BOY- ( 0  -  I  )*2 
I  OF  -  J 

GO  TO  6543 
ELSE 

GO  TO  5432 
END  IF 
ENO  IF 

100  FORMAT! 2! 3F 10.0, 21 1 > ) 

201  FORMAT! 2! 3F 10. 6 , I  1 ) > 
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6543  TYPE*. 'WALL  SECTIONS  ARE  THOSE  FOR  WHICH  WAKE  EFFECTS  ARE  ABSENT’ 

TYPE*. 'SAY  0  IF  THERE  ARE  NO  WALLS’ 

ACCEPT*, IWALL 
I F { IWALL . EQ . 0 )  GO  TO  9B76 
TYPE*. 'HOW  MANY  WALL  SECTIONS?' 

ACCEPT*, NWALL 

IF ( IWALL . EO. 1)  WALL*. TRUE. 

9876  IF(MOO(NONBOY,2) .EO.0)  THEN 

N  *NON8D Y/ 2 
EVEN* -TRUE . 

ELSE 

N  =  <  NONBDYM  >/ 2 
EVEN* . FALSE . 

END  IF 
C 

C  NOW  ARRANGE  THE  ONBOOY  DATA  ALONE  ONE  COORDINATE  PER  LINE. 

C 

REWIND  2 
DO  10  1*1 ,N 

WRITE  1  2)  X1(1).Y1(I), 21(1), LMI), Mill) 

IF( I .EQ.N. AND. EVEN. EO.  .FALSE.  )  GO  TO  10 
WR I TE ( 2 )  X2C I  )  ,Y2( I )  ,Z2( I  I  ,L2( I ) ,MZ< I ) 

10  CONTINUE 

LASSTR 1 = . FALSE . 

LASSEC*. FALSE . 

TYPE*, 'ARE  THERE  FUSELAGE  SECTION?  SAY  0  IF  NONE' 

ACCEPT* . I F  USE 

IF ( I  FUSE . EO .0)  FUSE*  .FALSE. 

IF< IFUSE.NE.0)  FUSE*  .TRUE. 

C 

C  NEXT  READ  THE  ARRANGED  DATA  TO  COMPLY  WITH  THE  SUBSEQUENT  PATTERN. 

C 

N*NONBDY 
REWIND  2 
DO  3210  1*1 ,N 

3210  READ  <  2  )  X( I ) , Y( I ) ,Z< I  )  ,L< I > ,M< I ) 

NEWSEC*0 
JJU-1 
I  TER  =0 
NSUN l =0 
KWALCT=0 
L I FT=0 

inn  CONTINUE 

K  1  *0 

DO  110  1  =  1  ,N 

IF (  I . LT. JJJ  )  GO  TO  1 10 

IF ( ITER . EQ. 3  !  GO  TO  112 

IF  1 L ( I > . EQ. 2 . AND .Ml  1 > .EO. 1 >  THEN 

GO  TO  111 

ELSE 

END  IF 

IF ( L (  I  ) . EQ. 2 )  SECT*. TRUE. 

IF (L< I ) .NE .2)  SECT*. FALSE. 

IF ( SECT )  NEWSEC  » NEWS EC  + 1 

IF (WALL . EQ .. TRUE .  .AND . NEWSEC . L E . NWAL L )  NSUN 1  - NSUN t  +  l 
IF ( NEWSEC . EO . 2 .AND . ITER . EQ. Z . AND . I . GT.JJJ )  GO  TO  111 
IF ( NEWSEC . EQ . 2 .AND . ITER . EO. 1 .AND . I . GT. JJJ )  GO  TO  2000 
IF ( NEWSEC .EQ. 2 .AND . ITER . EQ.0)  GO  TO  1000 
NEWSEC* 1 
I F ( SECT )  N 1 *0 

IF ( L (  I  )  . GE . 1 )  NEWSTR1-  .TRUE- 
IF(L<  I  ) .EQ.0)  NEWSTR 1  *  .FALSE. 

I F ( NEWSTR 1 )  GO  TO  1001 
GO  TO  1002 

1001  N 1 *N I ♦ 1 

IF ( NEWSEC .EQ. 1 )  NLINES1-N1 
GO  TO  1003 

1002  I F ( NEWSTR 1 )  MM1*0 
MM  1 *MM 1  +  1 
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IF < NEUSTR 1 ,E0. .FALSE. . AND . NEWSEC . EQ . 1 )  ML  I NES 1 -MM  1 
IF (SECT. EO.  .TRUE.  . AND . NEWSEC . EQ .  1 )  1 0-1 
K- I  -  I  0* 1 

IF<LASSTR1 .EQ. .FALSE. >  GO  TO  1006 
K  1  -<  1  +  1 

X  LAST  1  (  K  1  ) -X  (  I  ) 

YLAST 1 ( K  1  > - Y  <  I ) 

ZLAST1 { K I  ) -2  < I > 

IF < K1 .GT . 1 .OR. ITER .EQ. 1 )  GO  TO  1009 
GO  TO  1008 

1003  IF i NEUSTR 1 .EO. .TRUE . .AND .  L  {  I +  MM1 ) . EQ . 2 . AND . NL I NES 1 . GT . 1 > 

1  LASSTR 1- .TRUE . 

GO  TO  1002 

1006  IF(K.EQ.l)  GO  TO  1007 
GO  TO  1010 

1007  XLEl-X(I) 

YLEl-Yt I ) 

ZLE1»Z( I ) 

GO  TO  1010 

1008  BASEW*ABStZLASTl(Kl l-ZLEl ) 

1009  XLAST IP 1 < Kl) -XLAST 1 (  K 1 1+0. 3*6 . *BASEU 
XLAST 1  P  2  (  K  1 ) -XL AST  IP  1 ( Kl  ) +0 . 7*6 . * BASEW 
YLAST 1 P 1 ( Kl ) -YLAST 1 ( Kl ) -0. 3* ( YLAST 1 ( Kl (-YLE1 ) 

YLAST 1 P  2 ( K 1  1-YLE1 

ZLAST1P1 ( Kl 1-ZLAST1 (Kl  ) -0 . 3* ( ZL AST  1 ( K 1 ) -ZL E 1 ) 

ZL AST  1 P  2 ( K 1  l-ZLEl 

1010  CONTINUE 

XOL  D ( N 1 , MM  1 ) -X ( I) 

YOL  D ( N 1 , MM  1  >-Y( I) 

ZOL  D ( N 1 , MM  1  )-Z<  I  ) 

GO  TO  110 

112  LIFT-LIFT+1 

WRITE(5)X( I ) ,Y(I >,Z( I).L(I).M( I ) 

110  CONTINUE 

111  CONTINUE 
I  TER  =  I  TER  +  1 
LASSEC- . TRUE . 

I F ( ITER.GT.3)  GO  TO  4321 
C 

C  NOW  REARRANGE  THE  INPUTS 

C 

141  CONTINUE 

IFIWALL.EQ.  .TRUE.  . AND.NSUN1.LE. NWAL L )  NLINES1-NLINES1-2 

DO  101  N 1  -  1 , NL I NES 1+2 
NEWM 1B-NLINES1+3-N1 
NEWM1T-NL INES1+1+N1 

IF (WALL . EQ. . TRUE. .AND . NSUN1 . LE . NWALL )  NL I NES 1 -NL I NES 1 +2 
IF(FUSE)  GO  TO  149 
J  1  • ( ML  I NES1  +  1  1/2 

149  IFIFUSE.EQ. .TRUE. .AND. NEWSEC. LE. 2)  J1-MLINES1 

IF (  I  TER . EQ. 1 )  NL 1*0  1 

IF( ITER.E0.2)  NL 2-0  1 

IF( ITER.EQ.3)  NL3-01 

00  103  00-1 ,01 

IF(FUSE)  GO  TO  150 

ML  INI (JO ) - ( ML  I NES 1+3  >/2-00 

ML IN1P (00 ) -ML INES1+1-0J 

ML  INI (J 1  i-ML IN1P (01 ) 

GO  TO  151 

150  I F ( NEWSEC . EQ . 1 )  ML  I  N  1  P  (  JJ  > -0  J 

I F ( NEWSEC . EQ . 2 )  ML  I N 1 ( JO ) -ML  I NES 1 -0 J  + 1 

151  IF ( N 1  .EQ.NLINES1  +  1  )  GO  TO  1011 
IF ( N 1 .EQ.NLINES1 +2  )  GO  TO  1012 
IF(FUSE.EQ. .TRUE. .AND. NEWSEC. EQ. 1 )  GO  TO  152 
XNEWIBtOO , NEWM IB  >*XOLD( N1 .ML  INI ( JO > ) 

YNEW1B! JJ .NEWM1B )-YOLD(Nl .ML  INI (OJ )  ) 

Z  NEW1 B ( 0 J , NEWM 1B)-ZOLO(N1 .ML  INI  I  00  )  ) 

IF(FUSE)  GO  TO  102 

152  XNEW1TI JO , NEWM 1 T ) -XOL  0 ( N 1 .ML  I  NIP (00 ) ) 
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YNEW1TI00 .NEWM1T  >-YOLD( N1 .ML  I  NIP <00  >  > 

ZNEW1 T( JO , NEWM 1T)*ZOLO(N1 , ML  INI P (00! > 

GO  TO  102 

1011  IF(FUSE.EO. .TRUE. .AND. NEWSEC. EQ.l)  GO  TO  153 
XNEW1 B  {JO  .  NEWMIB  I-XLAST1P  1  (MUNI  <  00  >  ) 

YNEV1 B ( 00 , NEWMIB 1-YLAST1P 1 (ML  INI (00  > ) 

ZNEW1 B ( 00 , NEWM1B  1-ZLAST1P 1 (ML  INI  (00  )  ) 

IF(FUSE)  GO  TO  102 

153  XNEW1 T ( 00 , NEVM 1 T I »XLAST 1 P1(MLIN1P(00 >  > 

YNEW1T(00,NEWM1T)«YLAST1P 1 (ML  IN  IP (ODD 
ZNEV1 T ( 00 , NEWM1 T >«ZLAST!P1(MLIN1P(DJ>  ) 

GO  TO  102 

1012  IF (FUSE. EO. .TRUE. -AND. NEWSEC. EQ. 1!  GO  TO  IS4 
XNEV 18(00 , NEWM1B ) -XLAST 1 P 2 ( ML  I N 1 (001) 

YNEW 1  B  (  00  ,  NEUM  1  B  )  *YLAST  1  P  2  !  ML  I  N1  (00  )  ) 

ZNEW 1 B ( 00 . NEWM1B  J-ZLAST1 P2( ML  INI (00 )  ) 

I F ( FUSE )  GO  TO  102 

LS4  XNEW1T ( 00 . NEWM 1 T  > -XLAST 1 P 2  I  ML  I  NIP  <00 )) 

YNEW1 T ( 0 0 , NEWM 1T)*YLAST1P2(MLIN1P<00>  ) 
ZN£W1T(00,NEWM1T)-ZLAST1P2(MLIN1P(00>  ) 

102  CONTINUE 

103  CONTINUE 

101  CONTINUE 

IF < FUSE. EQ. .TRUE . .AND . NEWSEC .EO. I )  GO  TO  11111 
IF(FUSE.EO. .TRUE. .ANO. NEWSEC. EQ. 2)  FUSE*  .FALSE. 

GO  TO  121 

1000  NEWSEC*  1 

1  TER  *  I 

lasstri*. false. 

000*000*NLINES1*MLINES1 
REWIND  5 
GO  TO  141 

2000  NEWSEC *2 

LASSTRI=.FALS£. 

0.00»000  +  NLINES1*MLINES1 
I  TER* I  TER  *  1 
GO  TO  141 

121  CONTINUE 

DO  130  1-1,01 

IF(WALL.EO. .TRUE. . AND . NSUN 1 . L E . NWAL L )  NL I NES 1 -NL I NES1 -2 
DO  140  0*1 , ( 2*NL I NES 1 +3 ) 

IFIWALL.EQ.  .TRUE.  .ANO.NSUN1  .LE.NWALL)  NL  I  NES  I  *NL  I  NES  I  -  •' 
IF(J.EQ.l)  LNEV1B! I , Jl-l 
IF(0 .EO. 1 .AND. I .EQ. 1)  LNEW1B(1.0)*2 
IF(0.GT.NLIN£Sl+2)  GO  TO  160 

WRITE(5>XNEW1B! I , J ) , YNEW1 8 ( I .0 ) .ZNEW1BI I ,0  > .LNEW1BI I .0 ) 

1  .MNEWIB(I.O) 

GO  TO  140 

160  WR I TE ( 5 1XNEW1 T ( I ,0 ) ,YNEW1T( 1 .0  > ,ZNEW1T( 1,0), LNEW1T ( 1,0! 

I  , MNEW 1 T ( I . 0 ) 

140  CONTINUE 

130  CONTINUE 

C 

C  FOR  MULTIPLE  SECTIONS 

C 

I F ( IWALL.EQ. 1 .AND. NWAL L.EO. NSUN 1)  WALL-. FALSE. 
IFILASSEC.EQ. .FALSE. )  GO  TO  11111 
000-000*NLINES1*MLINES1 
I F  < ITER . EQ. 3  )  GO  TO  11111 

4321  IF ( IFUSE . EO. 1 . OR . NWALL .EO. 1 )  I WAKE *4* ( NL2+NL3) 

I F ( IFUSE.EQ.0.AND. IWALL .EO.0)  IWAKE-4* < NL 1 *NL2+NL 3 > 

I F  <  NWAL L . EQ . 2 )  IWAKE- 4 * { NL 3 ) 

NN*U0O  *■  IWAKE  +  L  I  FT 
I I-MOD! NN , 2  > 

I F ( 1 1 . Nc . 0  )  NN-NN+1 
C 

C  NEXT  REWRITE  TWO  PANELS  PER  LINE 

C 


REWIND  5 
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o  o  o 


40 

1 

194 
123 
2  00 
C 
C 
C 

1 


22222 


33333 


00  40  I-l.NN/2 

READ! 5)  XNEW1 ( I ) , YNEW1 ( I ) .ZNEW1 ( I ) , l NEW 1 ( I ) , MNEW1 ( 1 > 

IF < I  I .NE .0.AND . I  . EO. NN/2)  GO  TO  40 

READ  <  5 ) XNEW2 ( I ) . YNEW2! I  1 , ZNEV2  < I ) ,LNEV2< I ) ,MNEV2< I ) 
CONTINUE 

WR ITE! 3, 200)  <  XNEW1 < I ) , YNEW1 ( 1 ) .ZNEV1 ( I ) . LNEW1< I ) ,MNEW1< I ) 
XNEV2< I ) , YNEW2 ( 1 > ,ZNEW2< I > , LNEW2< I ) , MNEW2 ( I >  ,  1-1 
CONTINUE 
CONTINUE 

FORMAT ( 2 ( 3F 10.6,2111) 

WRITE  DOWN  THE  OFF  BODY  POINTS  IN  THE  SAME  OUTPUT  FILE. 

WRITE! 3.201)  (XI ( I ) ,Y1 ( I ) ,Z1 ( I ) ,L1 (I)  ,X2< I ) ,Y2l 1 ) ,Z2( I > ,L2 
I " I  ON  +  1 , IOF  ) 

FINALLY  READ  AND  WRITE  THE  BOUNDARY  LAYER  PARAMETERS. 

READ (  1 ,22222 )  UI ,R I , ISM. I  BETA , NOCAL , NTRANS 
WRITE (3. 22222)  UI  ,RI .  ISM. I  BETA. NOCAL .NTRANS 
FOR  MA  T(2F15.6,3I5, 12) 

READU, 33333)  (NXT(I).I-l.  NTRANS  > 

WRITE(3, 33333)  <  NXT< I > , I  *  1 .NTRANS ) 

FORMAT! 15 ) 

STOP 

END 


!  NN/2  ) 


(  1  )  . 
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C*OECK  NOLIFT 
C  LIST. NONE 

SUBROUTINE  NOLIFT 

.....THIS  SUBPROGRAM  COMPUTES  THE  GEOMETRIC  QUANTITIES  OF  THE 
*****NON-L  IFTING  ELEMENTS.  REFER  TO  REPORT  NO.  E.S.  40622. 


C 

C 

C 


REAL 

NX 

.NY  . NZ , 

IXX. 

IXY,  IYY 

DIMENSION 

XIN 

4  ) 

YIN  I  4  ) 

DIMENSION 

ZIN 

4  ) 

XA  I  500  ) 

DIMENSION 

YA  ( 

500  > 

ZA  (  500  > 

DIMENSION 

XB  < 

500  ) 

YB  (  500  ) 

DIMENSION 

ZB  ( 

500  ) 

XI  I  4  ) 

DIMENSION 

ETA 

4  ) 

TITLE  (  18  > 

DIMENSION 

XC  I 

1100  ) 

YC  (  1100  ) 

DIMENSION 

ZC  ( 

1100  ) 

XN  I  1100  ) 

DIMENSION 

YN  { 

1  100  ) 

ZN  (  1100  > 

DIMENSION 

UNITVX  (  1100) 

UNITVY  (  1100 

DIMENSION 

UNITVZ  (  1100) 

A  (  1100  1 

DIMENSION 

NSORNL ( 50) 

COMMON 

/ 

POINTS 

/ 

XA.  YA 

,  ZA, 

XB.  YB.  ZB 

COMMON 

/ 

CONFLG 

/ 

TITLE, 

CASE 

.  LIFSEC,  LIST. 

KUTTA, 

MPR 

COMMON 

/ 

MISC 

/ 

N.  LABEL.  NPRT,  LPRT.  MMIN 

COMMON 

/ 

IADDRS 

/ 

I  SAVE . 

OJ  . 

I  COUNT 

I  OUT .  KONTRL, 
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COMMON 
COMMON 
COMMON 
COMMON 
COMMON 
DATA  BTYPE 


/  CTABLE  /  XC.  YC,  ZC  .AMACH 
/  NORMAL  /  XN,  YN,  ZN 

/  ATABLE  /  A 

/  UVETOR  /  UNITVX,  UNITVY.  UNITVZ 

/  SCALE  /  FC.  BL 


/  4HNLIF  / 


C 

C 

CT500  FORMAT  ( 
CT  1 


CT  2 
CT  3 
CT  4 
CT  5 
CTS10 
CT  1 
CT520 
CT  1 
CT530 

500 

I 

501 

502 
C 

C 

c 


M,  7X.  4 (  1HX,  1 1 X  ).  ZHNX ,  10X.  SH  X0  . 
7HTYPE  OF  /  27X,  4!  1HY,  11X  ).  2HNY,  10X. 
1HT,  10X,  7HELEMENT  /  27X.  4<  1  HZ  .  1 1 X  ), 

10X,  1  HA  /  MX.  6H-- 
10<  1H-  >.  4X,  9  !  1H- 


)  , 


.  4 !  3X. 

4X.  10<  1H- 


1H0,  14X.  5HN 
10X,  1  HD .  10X 
5H  Y0  ,  10X, 

2HNZ ,  10X.  5H  Z0 

9H -  ),  3X  . 

.  5X.  7 {  1H-  )  /  ) 

1H0,  15X,  M.  4F12.6,  2F13.S,  1PE14M  / 

20X ,  0P4F12.6,  2F13.6.  1PE14M  >  ) 

1H0,  1  I X .  214,  4F12.6,  2F13.6,  1  PE  14  M ,  6X .  A4  /  ( 

0P4F12.6,  2F13.6,  1PEM.4  )  ) 

1H0.  S3X.  13UH-I,  3H  *  .  14I1H-1  ) 

FORMAT { 12X , 4HN  M.5X , 2HX0, 8X , 2HY0 , 8X . 2HZ0 , 8X . 2HXN , 8X , 2HYN , 
8X,2HZN,9X,1H0.9X.1HT,9X,1HA,6X,4HTYPE> 

FORMAT! 10X . 2 13 . 9F 10. 6 , 3X ,A6 ) 

FORMAT!  1  3X . 13. 9F 10. 6  ) 


FORMAT  ! 

( 

FORMAT  ( 


FORMAT  ! 


20X  . 


SUNBET-SORT! 1 . -AMACH'AMACH ) 
WRITE  (  4  )  LABEL 
DO  2000  I  -  1.  MM  IN 


KONTRL  - 

KONTRL  ♦  1 

XIN 

!  1  )  - 

XA 

!  I  ) 

XIN 

(2)  - 

XB 

(  I  > 

XIN 

(3)  - 

XB 

!  1*1  ) 

XIN 

(4)  - 

XA 

(  I  ♦  1  > 

i 

YIN 

!  1  )  • 

YA 

1  I  ) 

YIN 

I  2  >  = 

YB 

(  I  ) 

YIN 

(3)  • 

YB 

(  1*1  ) 

\ 

YIN 

(  4  )  - 

YA 

(  1  +  1  ) 

ZIN 

(  1  )  ■ 

ZA 

(  I  ) 

ZIN 

(2)  • 

ZB 

(  I  ) 

ZIN 

(  3  )  * 

ZB 

(  1*1  ) 

i 

ZIN 

!  4  )  - 

ZA 

(1*1) 

.  DO  1 

500  K 

-  1 

,  4 

1  XIN 

(  K  ) 

m 

XIN  !  K 

) 

/ 

FC  f 

|  YIN 

!  <  ) 

*  ' 

YIN  !  K 

> 

/ 

FC  I 

1  ZIN 

(  K  ) 

» 

ZIN  (  K 

) 

r 

FC  F 

600  CONTINUE 


•FORM  DIAGONAL  VECTORS.  EQUATION  (  64  ) . 

T1X  -  XIN  (3)  -  XIN  (1) 

T2X  -  XIN  (4)  -  XIN  ! 2 ) 

T1Y  -  YIN  ! 3 )  -  YIN  ! 1 ) 

T2Y  -  YIN  ! 4 )  -  YIN  ! 2  > 

T1Z  -  ZIN  (3)  -  ZIN  ( 1 ) 

T2Z  -  ZIN  (4)  -  ZIN  (2) 

•FORM  CROSS  PRODUCT  N  •  TZ  X  Tl.  EQUATION  !  65  > . 

NX  -  T2Y  *  T 1 Z  -  T1Y  *  T2Z 

NY  -  T1X  •  T2Z  -  T2X  *  T1Z 

NZ  -  T2X  •  T1Y  -  T1X  *  T2Y 

VN  •  SORT  !  NX  *  NX  *  NY  •  NY  ♦  NZ  •  NZ  ) 

•FORM  UNIT  NORMAL  VECTOR,  EQUATION  I  66  >. 
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NX  -  NX  /  VN 
NY  •  NY  /  VN 
NZ  -  NZ  /  VN 
C 

c.*.. -COMPUTE  AVERAGE  POINT.  EQUATION  (  68  )  . 

C 

AVX  -  0.25  *  (  X  I N  <  1  )  ♦  X  I N ( 2  >  ♦  X  I N  <  3  >  ♦  X I N  <  4  >  ) 

AVY  •  0.25  *  <  YIN(l)  ♦  YIN<2>  ♦  VINO)  +  Y2N(4)  ) 

A VZ  -  0.25  *  (  ZIN<1)  ♦  ZIN<2>  ♦  ZIN(3>  ♦  Z I N  <  4  >  > 

C 

C  *  *  *  *  *COMP  UTE  PROJECTION  DISTANCE.  EQUATIONS  (  69  )  AND  (  71  ). 

C 

0  -  NX  *  ( AVX-X I N  <  X  > }  +  NY  *  (AVY-YIN(l))  ♦  NZ  *  (AVZ-ZIN(l)) 
P  D  «  ABS  I  0 ) 

C 

C*.««. EQUATIONS  (73)  AND  (74). 

C 

T  •  SORT  (  T1X  *  T1X  +  T1Y  *  T1Y  ♦  T1Z  *  T1Z  > 

T1X  «  T1X  /  T 
T1Y  =  T1Y  /  T 
T1Z  *  T1Z  /  T 
C 

C*****EQUATION  (  75  >. 

C 


T2X  * 

NY 

*  T1Z 

NZ 

*  T1Y 

T2Y  - 

NZ 

*  T1X 

NX 

*  T1Z 

T2Z  * 

NX 

*  T1Y 

NY 

*  T1X 

C 

C-«...C0MPUTE  COORDINATES  OF  CORNER  POINTS  IN  REFERENCE  COORD.  SYSTEM. 
C*«*«*EQUATION  (  72  ) 

C 

C 

C  -  « - 

C 

DO  1000  J  -  1 .  4 

XP  -  XIN  (J)  +  NX  *  D 

YP  -  YIN  (0)  ♦  NY  *  D 

ZP  •  ZIN  (J)  +  NZ  *  D 

D  »  -0 

XDIF  -  XP  -  AVX 
YD  IF  -  YP  -  AVY 
ZD  I F  »  ZP  -  AVZ 
C 

C«****TRANSFORM  CORNER  POINTS  TO  ELEMENT  COORDINATE  SYSTEM  (  XI,  ETA  ) 
1  C-----VITH  AVERAGE  AS  ORIGIN.  EQUATION  (  80  ) . 

I  C 

1  X I ( 0  )  -  T1X  *  XDIF  +  T1Y  *  YD  I F  *  T1Z  *  ZDIF 

1000  ETA(0>-  T2X  *  XOIF  +  T2Y  *  YD  I F  ♦  T2Z  »  ZDIF 

C 

C  - * - 

C 

C 

C..*»«C0MPUTE  CENTROID.  EQUATION  (  81  ). 

1  C 

XI0  -  0.3333333E0  •  (  XII4I  *  (  ETA(l)  -  ETA( 2 )  )  ♦  XI  (2) 

1  *  (  ETA ( 4 )  -  ETA( 1 >  )  )  /  (  ETA ( 2 )  -  ETA(4)  ) 

ETA0--0 . 3333333E0  *  ETA(l) 

C 

C«***»OBTAIN  CORNER  POINTS  IN  SYSTEM  WITH  CENTROID  AS  ORIGIN,  EQ.  (  82 
C 

DO  1020  J  •  1 ,  4 
XI  (0)  *XI  (0)  -XI0 
1020  ETA( J  >  »  ETA(J)  -  ETA0 
C 

C»««. -COMPUTATION  AIDS. 

C 

ETAZM1  •  ETA  (2)  -  ETA  (1) 

~  ETA3M2  *  ETA  (3)  -  ETA  (2) 

ETA4M3  *  ETA  (4)  -  ETA  (3) 
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oooo  ooo  ooo  ooo  ooo  non 


ETA1M4  *  ETA  (1)  -  ETA  (4) 

XI  1M2  -  XI  < 1 )  -  XI  (2) 

XI2M3  •  XI  (2)  -  XI  (3) 

XI3M4  -  XI  (3)  -  XI  (4) 

XI4M1  -  XI  (4)-  XI  ( 1 > 

ETA2P4  -  ETA  (2)  ♦  ETA  (4) 

XI3M!  -  XI  { 3 )  -  XI  (1) 

XI4M2  -  XI  (4)  -  XI  (2) 

ETA2M4  «  ETA  (2)  -  ETA  (4) 

XI 1234  »  XI  (1)  ♦  XI  (2)  +  XI  (3)  ♦  XI  (4) 

C 

C****« TRANSFORM  CENTROID  TO  REFERENCE  COORDINATE  SYSTEM.  EQ.  (  83  >. 

C 

C 

C  -  * - 

C 

XCENT  -  AVX  *  T1X  *  XI0  ♦  T2X  »  ETA0 

YCENT  -  AVY  +  T1Y  *  XI0  +  T2Y  *  ETA0 

ZCENT  -  AVZ  ♦  T1Z  *  XI0  +  T2Z  *  ETA0 

•COMPUTE  LARGER  DIAGONAL  VECTOR,  EQUATION  t  84  ) . 

TSQ  -  AMAX1  (  XI3MI  •*  2.  X 1 4M2  *•  2  ♦  ETA2M4  ••  2  ) 

T  -  SORT  (  TSO  ) 

•COMPUTE  AREA,  EQUATION  (  85  > . 

AREA  -  0.5  *  XI3M1  •  ETA2M4 

•COMPUTE  SECOND  MOMENTS  IXX.  IXY,  IYY.  EQS .  (  86  >  TO  (  88  >. 

IXX  -  8.333333E-2  *  XI3M1  *  (  ETA  (1)  •  XI4M2  •  XU234  +  ETA2M4 

1  *  (  XI  (1)  *  (  XI  (1)  +  XI  <  3  >  )  ♦  XI  (3)  *•  2  )  ♦  XI  (2) 

2  •  ETA  (2)  •  (  X I  1 234  -  XI  (4)  )  -  XI  (4)  *  ETA  (4) 

3  *  (XI 1234  -  XI  (2)  )  ) 

IXY  *  4.166667E-2  •  XI3M1  *  <  2.8  •  XI  (4)  •  (  ETA  (1)  •«  2  - 

1  ETA  <41  •*  2  )  -  2.8  *  XI  (2)  *  (  ETA  (1>  *•  2  -  ETA  (2>*»2 

2  *  <  XIII)  ♦  XI C  3  >  >  *  ETAZM4  •  (  2.0  •  ETA  (1)  ♦  ETA2P4  >  ) 

IYY  -  8.333333E-2  *  XI3M1  *  ETA2M4  *  (  (ETA  (I)  ♦  ETA2P4 )  *•  2 

1  -  ETA  (I)  *  ETA2P4  -  ETA  (2)  •  ETA  (4)  ) 

•COMPUTE  CONSTANTS  FOR  EQUATIONS  (  42  )  AND  (  43  ) .  ALSO  EQ.  (  45 

0I2SQ  «  XI1M2  *•  2  ♦  ETA2M1  **  2 
D 1 2  -  SORT  (  DI2SQ  ) 

D23SQ  »  XI2M3  **  2  ♦  ETA3M2  ••  2 
D23  *  SORT  (  D23SQ  ) 

D34SQ  »  XI3M4  **  2  ♦  ETA4M3  **  2 
034  -  SORT  (  D34SO  > 

D41SQ  -  XI4MI  ••  2  ♦  ETA1M4  •*  2 
0 4 1  -  SORT  (  D41SQ  ) 


XC  <  KONTRL  >  »  XCENT 
YC  (  KONTRL  )  *  YCENT 
ZC  (  KONTRL  )  *  ZCENT 
UNITVX  <  KONTRL  )  •  T1X 
UNITVY  (  KONTRL  )  •  T1Y 
UNITVZ  <  KONTRL  )  •  TiZ 
OO  1500  K  *  1 ,  4 

XIN  (  K  )  *  XIN  (  K  )  •  FC 

YIN  (  K  )  -  YIN  (  K  )  •  FC 

ZIN  (  K  )  -  ZIN  (  K  )  *  FC 

1500  CONTINUE 


••••SPRINT  RESULTS  --  SECTION  9.4  THE  FIRST  OUTPUT. 


119 


CS  NOTE  THAT  THE  COMPR ESS  I B I L I  TV  CORRECTION  IS  BE  INS  MADE  HERE. 


XO 

XCENT 

•  FC 

YO 

YCENT 

*  FC 

/SUNBET 

ZO 

ZCENT 

*  FC 

/SUNSET 

TT 

T 

*  FC 

/SUNBET 

DD 

PD 

*  FC 

/SUNBET 

AR 

AREA  * 

'  F, 

*  FC  /<SUN8ET*SUNBET) 

IF  < 

NPRT  .GE. 

1  1  ) 

GO  TO  1750 

IF (NPRT.GT . 38 >  GO  TO  1750 

NPRT  -  NPRT  +  1 

IF  (  I  .EQ.  1  )  GO  TO  1760 

CT  WRITE  (  6,  510  )  I.  XIN.  NX.  XO ,  DD .  YIN.  NY,  YO.  TT .  ZIN, 

CT  1  NZ,  ZO,  AR 

WRITE  (  6,  502  1  I .XO.YO.ZO.NX.NY.NZ.DO.TT.AR 
GO  TO  1770 
1750  NPRT  -  0 

CALL  HEADER 

CT  WRITE  !  6.  530  ) 

WRITE  <  6,  500  ) 

CT760  WRITE  (  6,  520  )  N,  I,  XIN.  NX.  XO,  DD.  8TYPE.  YIN.  NY,  YO, 

CT  1  TT,  Z1N,  NZ,  ZO.  AR 

1760  WR ITE  f  6 , 501 )  N , I , XO , YO , ZO . NX . NY , NZ , DD , TT , AR . BTYPE 
1770  XN  (  KONTRL  )  -  NX 
YN  <  KONTRL  >  «  NY 
ZN  I  KONTRL  )  *  NZ 
A  I  KONTRL  )  *  AREA 
C 

c«*.«  WRITE  28  QUANTITIES  ON  TAPE  4  AS  ONE  LOGICAL  RECORD. 

C 

2000  WRITEU)  XCENT ,  YCENT ,  ZCENT ,  T1X,  T1Y,  T1Z,  T2X ,  T2Y ,  T2Z ,  NX, 
1  NY,  NZ,  XI(l),  ETA(l),  X  I ( 2 ) ,  ETA (  2  I  ,  XI(3).  ETA( 3  > , 

Z  X I <  4 ) .  ETA<  4 ) ,  TSQ .  AREA.  IXX,  IXY.  IYY,  D12.  D23. 

3  034,  D4I 

I COUNT  •  ICOUNT  +  MM IN 

RETURN 

END 


120 


SUBROUTINE  BSETUP 


C 

c 

c 

c 

c 

c 

c 

c 


c 


c 


THIS  SUBROUTINE  SETS  UP  ALL  THE  NECESSARY  PARAMETERS  AND 
REQUIREMENTS  BEFORE  CALLING  THE  BOUNDARY  LAYER  PROGRAM.  ALSO 
COMPUTES  THE  NEW  SIGMA'S,  THEN  LOOPS  BACK  TO  SUBROUTINE  COMFLO 
TO  COMPUTE  THE  FINAL  SOLUTIONS. 


NOTE  HERE  THAT  THE  NECESSARY  CHANGES  FOR  THE  INTRODUCTION 
OF  VISCOUS  EFFECTS  ON  THE  NON-LIFTING  BODY  ALSO  IS 
INCLUDED.  ALSO,  SOME  INPUT  MODIFICATIONS  ARE  MADE. 


DIMENSION 

BLINDI2) 

DIMENSION 

NTYPE 

I  50 

)  , 

NLINE  i 

[  50  ) 

DIMENSION 

TITLE 

(  18 

>  , 

NSORCE  i 

[  10  )  , NSORNL ( 30 

DIMENSION 

NWAKE 

<  10 

)  , 

NSTRIP  i 

I  10  ) 

DIMENSION 

IGN  ( 

50.10 

)  , 

1G1  (  50,10  ) 

DIMENSION 

KOUNT 

(  50 

)  , 

NLT 

:  500  ) 

DIMENSION 

UNITVX 

(  1  100 

)  , 

UNITVY 

!  1100  ) 

DIMENSION 

UNITVZ 

(  1  100 

>  , 

VELX 

I  1100  ) 

DIMENSION 

VELY 

(  1100 

)  , 

VELZ 

1100  ) 

DIMENSION 

XCENT 

(  1  100 

)  , 

YCENT 

1  100  ) 

DIMENSION 

ZCENT 

(  1  100 

)  , 

XC 

100  ) 

DIMENSION 

YC 

(  100 

)  , 

ZC 

100  ) 

DIMENSION 

X8 

<  100 

)  . 

YB 

100  > 

DIMENSION 

ZB 

<  100 

)  , 

V 

100  ) 

DIMENSION 

DEL 

(  100 

>  , 

SM 

100  > 

DIMENSION 

XINP 

(  500 

)  , 

YINP 

500  ) 

DIMENSION 

ZINP 

<  500 

>  , 

XLINP 

50  ) 

DIMENSION 

YL  INP 

(  50 

)  , 

ZLINP 

50  > 

DIMENSION 

RHS 

(  1  100 

>  . 

DQDS 

100  ) 

DIMENSION 

VEL 

(  1  100 

)  , 

DOS 

100  ) 

DIMENSION 

VB 

(  100 

)  , 

D  (  100  ) 

DIMENSION 

T 

(  100 

)  , 

BT ( 100) , 

S  (  100  ) 

DIMENSION 

DELS 

(  100 

)  , 

IXFLAG 

10  ) 

DIMENSION 

PCHORD 

(  51 

>  , 

SK 

51  ) 

DIMENSION 

VK 

(  51 

)  9 

XK 

51  > 

DIMENSION 

ZK 

(  51 

)  , 

DELK 

51  ) 

DIMENSION 

DF 

(  51 

)  , 

TF 

1  51  ) 

DIMENSION 

BETAK 

(  51 

) 

COMMON  / 

CONFLG 

/ 

TITLE, 

CASE,  LIFSEC,  LIST,  NOFF,  IOUT,  KONTRL . 

KUTTA, 

MPR 

COMMON  / 

INFORM 

/ 

NSORCE 

,  NWAKE, 

NSTRIP.  IG1. 

IGN 

COMMON  / 

CTABLE 

/ 

XCENT , 

YCENT , 

ZCENT  .AMACH 

COMMON  / 

IMPORT 

/ 

NON,  NFLOW 

COMMON  / 

BO  INFO 

/ 

KOUNT , 

NLINE , 

NTYPE.  NLT ,  ISECT 

COMMON  / 

VBLOCK 

/ 

VELX  , 

VELY,  VELZ 

.  VEL 

COMMON  / 

UVETOR 

/ 

UNITVX 

,  UNITVY, 

UNITVZ 

COMMON  / 

RSPARM 

/ 

12,  SM 

,  VB,  DEL, 

DOS 

COMMON  / 

NL  INP 

/ 

NL.  XINP,  YINP, 

ZINP 

COMMON  / 

LTINP 

/ 

LT,  XLINP,  YLINP 

,  ZLINP 

COMMON/BL 12/TI ,RMI 

,UI  ,RI ,PR .PRT.FK, 

RL.RMUI , R  HOI ,PSI .HE 

1  ,UE(  100  >  ,R0U00)  ,TV<  100  >  ,QW(  100)  ,RP<  100)  .  FW (  100  ) 

2  , BR< 100) ,TE< 100) ,RHOE< 100) ,RMUE< 100) ,GV< 100) ,GPU( 100) 

3  ,RF1  ( 100)  .RFZU00)  ,  YS  (  100 )  ,  IGX 1  (  100 )  .  FPW(  100  > ,  ROL  (  100 
COMMON  /  EXTRA  !  IXFLAG,  SKIP,  METHOD 

COMMON  /  APPA  /  NKK.NFF .MAT, 1ERR.KPRINT 
COMMON  /  NLFVIS  /  NLFSTR  ,  NSORNL  .  NONLIF 


DATA 


ZERO  /  0.0E0  I 


DATA  BLINO/2HLO.ZHUP/ 

DATA  I  ONE  /  1 

DATA  IRSIDE  /  16 

DATA  P CHORD  / 

1  0.011730, 

2  0.032724,  0.038007, 

3  0.077742,  0.089022, 


/ 

/ 

0.000000,  0.002273, 
0.014384,  0.017271 , 


0.044048, 
0. 101623, 


0.004555, 

0.020464, 

0.050945, 

0.115622, 


0.006867, 

0.024046, 

0.058796, 

0.131092, 


0.009244 , 
0.028102, 
0.067697 . 
0.148098, 
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nnnn 


c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


4 

5 

6 

7 

8 


0. 166692. 
0.313123. 
0.515000. 
0. 745470, 
0.943639 . 


0. 186918, 
0.343287. 
0.552599, 
0.782969 , 
0.966931  . 


0.208807, 
0.37498  1  , 
0.590864 , 
0.8  19281  , 
0.985936, 


0.232376  . 
0.4081 19 . 
0.629560, 
0. 853981  . 
1 .000000 


0.257630, 
0 .442592 , 
0. 66S427 , 
0.88661 1 , 

/ 


0.284555 

0.478271 

0.707171 

0.916674 


SUNBET-SQRT<  1 . -AMACH*AMACH  > 

10  FORMAT  <  1H0,  35X,  ' 1 NTERMED IATE  PRINT  IN  SUBROUTINE  SETUP’ 

( DIMENSIONS  IN  FEET) •  ) 


20  FORMAT 
1 

30  FORMAT 
1 
2 

32  FORMAT 

1 

2 

40  FORMAT 
99  FORMAT 

41  FORMAT 

50  FORMAT 

51  FORMAT 
1 

2 

52  FORMAT 

53  FORMAT 
1 

54  FORMAT 
1 

55  FORMAT 
60  FORMAT 
70  FORMAT 

1 

2 

80  FORMAT 
92  FORMAT 
1 

94  FORMAT 
1 

1 

95  FORMAT 

96  FORMAT 

97  FORMAT 


( 


1 H0 ,  4 1 X ,  ’SECTION  NO. 
•TOTAL  STRIPS  »  ’  .  13 


12.  4  X ,  ’TYPE 


1H0.16X.  1 BL .  STRIP 
’Y’ .  I5X,  ’Z’ .  15X. 
9 <  1 H- > .  3X  ,  8 (  1 H -  )  . 
1H0,  ax.  'STRIP  NO. 
'Z' .  15X.  ’V’ .  15X. 
9  <  1 H - > .  3X ,  8( 1H- ) . 
1H  .  1 9X ,  A5 .  10X. 

A5.  10X, 
12,  10X. 


■  ’ .  II.  4X. 

'  X  ’  .  1 5X  , 

/  1  7X  . 

’X’ ,  15X. 

’  DQDS ’  /  9X . 


’,  3X ,  'C.  POINT’.  9X , 

•V  ,  15X.  ’S' 

5<  3X ,  131 1 H- )  )  ) 

•  ,  3X ,  ’C.  POINT’  .  9X. 

•S',  14X.  ’DEL’.  13X. 

6  <  3X .  13UH-)  >  ) 

12.  6X .  E13.6,  4 (  3X .  E13.6  >  ) 

1H  ,  16X.  A5 ,  10X,  12.  6X ,  E13.6.  4(  3X.  E13.6  >  > 

1H  .  11X,  12,  10X.  12.  3X.6F16.6) 

1 H0 ,  53X.  ’THE  NEW  RIGHT  -  HAND  SIDE’  /  ) 

1H0,  24X.  1 HK ,  10X,  6HPCHORD.  10X ,  8HXK  (  K  ).  9X . 

8HZK  (  K  ).  9X ,  8HSX  (  X  ).  9X .  8  HVK  ( 


3H - .  9X  6 ( 1 H- ) ,  10X,  8 ( 1 H- > . 


1H  ,  22X, 
1H0,  37X. 
10HDEL  ( 
1H0,  37X, 
10HBETAK ( 
1H  .  36X. 
1H  ,  7X , 
1H0,  14X, 

14,  3X , 
’Z0  -  ’  . 
1H0,  46X 
1H0,  46X, 


13,  S<  3X ,  F14.6 
1  H  I  , 

I  ) 

1  HK  , 

K  ) 

13. 

8E15. 


)  ) 


K  )  /  24X, 
31  9X,  8< 1H-  >  )  / 


1 1 X .  8HVB  t 
5X .  141 1 H- > 
1 1 X .  8HVK  ( 
5X  ,  14  1 1 H- ) 


I  )  , 

>  / 

K  )  , 
)  / 


10X  , 


10X  , 


•STAGNATION  PANEL 
•  ,  El  3. 6.  3X , 


9X .  8HSM  <  I  ) . 

/  37X,  3H - .  3 ( 

9X ,  3HSX  <  K  >. 

/  37X ,  3H - ,  3 ( 

3 (  5X ,  E14.6  >  ) 

.6  ) 

,  'STRIP  NO.  =  ' .  12,  3X . 

'X0  *  ’ ,  E13.6.  3X ,  'Y0  - 
E13.6  ) 

,  'NO  SEPARATION  CAN  BE  FOUND  FOR  THIS  STRIP') 
•TOTAL  POINT  IN  SETUP  NOT  EQUAL  TO  TOTAL', 
’CONTROL  POINT.  PROGRAM  ENDS  FOR  CORRECTION' 
1H0,  4 7X ,  'ELEMENTS  INPUT  TO  B.L.  PROGRAM  =  ',  15  /46X 

'ELEMENTS  RETURN  TO  CALLING  PROGRAM  *  ’,  15  /  61X, 
'PROGRAM  STOPS'  /  ) 

2F15.0.  315  ) 

1H0,  4 7X ,  37HBOUNDARY  LAYER  PROGRAM  OUTPUT  FOLLOWS  /  1 

15  ) 


FIRST,  SET  UP  SOME  COUNTERS 

LS  --  LIFTING  SECTION  COUNT. 

NOS  --  NON  LIFTING  SECTION  COUNT. 

KK  --  ELEMENT  COUNT 

LI  --  INITIAL  STRIP  COUNT. 

L 2  --  FINAL  STRIP  COUNT 

KR  --  ORDER  FOR  ARRANGING  THE  S’S. 

KN  --  OVERALL  STRIP  COUNT 


Q*S,  DO  /  DS"S. 


I  ERR 
IM 

KSTR IP 
LS 
KK 
LI 

KTOTAL 

KN 

NL 


51 

0 

0 

0 

1 

0 

0 

0 


READ  <  5,  95  )  UI.  RI,  ISM.  IBETA,  NOCAL 


UI  -  REFERENCED  VELOCITY  (  FT  /  SEC  ) 

RI  ---  REYNOLD’S  NUMBER  PER  FOOT. 

ISM  --  0,  NO  SMOOTHING. 
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i 


C  IBETA  --  0,  NO  SPECIAL  BETA  CALCULATION. 

C  NOCAL  --  NON-ZERO.  NO  DELTA  STAR  GENERATED  FOR  NON-LIFTING  SEC. 

C 

C 

C  TAKE  ONE  SECTION  AT  A  TIME,  PROCEED  TO  PREPARE  FORMING  OF  A 

C  BOUNDARY  LAYER  STRIP. 

C 

DO  1 500  IS*  1  .  ISECT 
JTYPE  -  NTYPE  f  IS  ) 

JSTRIP  *  NLINE  (  IS  ) 

L  2  »  L  1  +  JSTRIP  -  1 

IF  (  JTYPE  .EQ.  1  )  GO  TO  500 
C 

C  JTYPE  IS  NON-LIFTING 

C 

KSTART  *  0 
K  1  *L  1 
LP-0 

C  LS»LS+1 

NS*NLF$TR 

C  IPOINT-NSORNL(NOS) 

REWIND  13 

DO  300  LK  *  LI ,  L 2 

NL  ■  NL  *  2 

KN  -  KN  ♦  1 

C  KSORCE  =■  IASS  (  NLT  (  KN  )  ) 

IF( IS.EQ. I .OR .LK.NE .L 1 )  KN-KN 
IF! IS.NE. 1 .AND.LK.EQ.L1 )  KN-KN+1 
L  L  * JSTR I P 

KSORCE  =  NSORNL I KN ) 

C  TYPE*. 'KSORCE, JSUNO ' . KSORCE , JSUNO 

J1  =■  KK  ♦  1 

J  2  =  KK  +  KSORCE 

IFIIS.LE.3)  NOCAL  =  1 
IFIIS.GT.3)  NOCAL  30 
IF  (  NOCAL  .GT.  0  }  GO  TO  240 
READ  I  5,  97  >  NXT 

12  «  KSORCE  ♦  1 

VB  ( 1 )  *  0. 

DELS! 1  >»  0. 

SM  <!>■  DELS  I  1  ) 

C 

C . THE  FOLLOWING  MODIFICATIONS  ARE  DONE  TO  INCLUDE  VISCOUS  EFFECTS 

C  ON  THE  NON-LIFTING  BODY  ALSO. 

DO  11111  KKK-1, KSORCE 
KK*KK+ 1 

B 1 l-UNITVXIKK) 

B12*UNITVY(KK> 

B13-UNITVZIKK) 

BVX*VELX( KK) 

BVY-VELYI KK ) 

BVZ*VELZ(KK) 

V(KKK)*VEL(KK) 

XCIKKKl-XCENTIKK) 

YC(KKK)*YCENT(KK) 

ZC(KKK)-ZCENTIKK) 

BT(KKK)*B11*BVX*B12*BVY+BVZ*B13 

11111  CONTINUE 
KMM 1  *  KSORCE- 1 

DO  11112  KKK  *  1 . KMM I 
BTPROD-BT(KKK)*8T(KKK+I ) 

IFIBTPROD.GE.ZERO>  GO  TO  11112 
KBSTRT-KKK 

IFtKBSTRT.LT. (KSORCE/2) 1LLLL-1 
IFIKBSTRT.EQ. ( KSORCE / 2 ) ILLLL-0 
1FILLLL.E0.1 >KBSTRT-KS0RCE/2 
GO  TO  11114 

11112  CONTINUE 
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c 


- H  0  STAGNATION  POINT  DETERMINED 

WRITE (6,80) 

GO  TO  4000 
11114  CONTINUE 

BDENOM=BT(KBSTRT*l ) -BT ( K8STRT ) 

IF(LLLL .EO.0)  GO  TO  12458 
X  I “XC ( KBSTRT ) 

Y I =YC < KBSTRT ) 

ZI=ZC< KBSTRT) 

GO  TO  13457 

124  58  XI  -  ( BT< KBSTR T»1  ) *XC ( KBSTR T ) - BT < KBSTRT > *XC < KBSTRT* 1 ) l/BDENOM 

YI*  ( BT ( KBSTRT  + 1 ) *YC ( KBSTR T )- BT ( KBSTR T ) *YC ( KBS TR T  + I )  )/ BDENOM 
ZI»  ( BT ( KBSTRT* 1 ) *ZC < KBSTRT 1  - BT ( KBSTR T ) *ZC ( KBSTRT* 1  )) /BDENOM 
13457  CONTINUE 

DO  11115  I8BL-1  .2 
CALL  HEADER 
WR ITE< 6 . 10) 

WRITEI6.20)  I S , JTYP  E , JSTR I P 
Y I »YI /SUNBET 
Z I *Z I /SUNBET 

WRITE(6,70)  LK, KBSTRT, XI ,YI ,Z1 

YI=YI*SUNBET 

Z I =Z I “SUNBET 

REWIND  13 

I8B“IBBL 

XB( 1 ) -X  I 

YB< 1 )“VI 

ZB<1 >-ZI 

IF ( IBB .EO. 2 )  GO  TO  11116 

I2=KBSTPT*1 

TYPE*,  ’  12’  ,12 

KM*  1 

DO  11117  K  =  2. 12 
KL  *  I  2 -KM 
XB(K)*XC(KL  ) 

YB( K  >»YC< KL ) 

ZB ( K ) *ZC ( KL ) 

VB< K )-V( KL ) 

KM“KM* 1 

11117  CONTINUE 
GO  TO  11118 
11116  CONTINUE 

I 2*KSORCE-KBSTRT* 1 
TYPE*. ’ 12  FOR  I BB*2 ’ , 12 
KM-KBSTRT 
DO  11119  K  =  2 , I  2 
KL  *KM* 1 


X B  <  K  >  -XC  (  KL  ) 

YB  (  K  >  -  YC  <  KL  > 

ZB ( K ) »ZC ( KL ) 

KM*KM* 1 

11119 

CONTINUE 

11118 

DO  I 1 109  KJ-2. 12 

DELX  *  <  XB  <  KO  )  -  XB  ( 

KO  -  I 

)  ) 

«*  2 

DELY  *  (  YB  (  KJ  )  -  YB  ( 

KO  -  1 

>  > 

**  2 

DELZ  *  (  ZB  (  KO  )  -  ZB  ( 

KO-t 

)  ) 

**  2 

DELS  (KJ)  -  SORT  (  DELX  ♦  DELY  ♦ 

DELZ  > 

SM  (KJ)  -  SM  (  KJ - 1  )  ♦  DELS  ( 

KJ  ) 

1  1109 

CONTINUE 

C 

XB  ( 1 )  *  0.5  *  (  XINP  (  NL 

)  ♦ 

XINP 

(  NL  - 1 

)  > 

C 

YB  (  1  )  -  0.5  *  (  YI NP  1  NL 

)  ♦ 

YI  NP 

(  NL- 1 

)  > 

c 

ZB  (  1  )  -  0.5  *  (  Z I  NP  (  NL 

>  ♦ 

Z  I  NP 

(  NL- 1 

)  ) 

c 

WRITE  (  6.  10  ) 

c 

WRITE  (  6,  20  )  IS.  JTYP  E  , 

JSTR I P 

c 

XB<  1  >*XB<  1  ) 

c 

YB (  I  )-Y8< I  )  /SUNBET 

c 

ZB (  1  >-ZB (  1  )  /SUNBET 

c 

WRITE  (  6,  70  )  LK,  KSTART 

,  XB 

(1  )  , 

YB  (  1  )  , 

ZB 

c 

XB( 1 > “XB ( 1 ) 

o  o  o  o  o  iinnn 


YB! 1 >-VB( 1 >*SUNBET 
ZB  < 1 >3ZSI 1  )*SUNB£T 
DO  100  K  J  3  2,  12 

_  k'k'  +  1 

XB  I  K  J  )  3  XCENT  (  KK  ) 

YB  ( K J ) 3  YCENT  <  KK  ) 

ZB  <  K J  >  =  ZCENT  (  KK  ) 

VB  (  K J  ) 3  VEL  (  KK  ) 

100  CONTINUE 

WRITE  I  S.  30  ) 

DO  102  I  3  1.  12 
XB I  I )=XB<  I  ) 

YB (  I ) = YB  <  I ) /SUNBET 
ZB ( I ) =Z8 ( I ) /SUNSET 
SMI  I  )  3SM I  I  ) 

WRITE  I  6.  40  )  BL TNDI IBBL > , I ,XB  I  I  )  .  YB  (n.ZBII),  VB  (I),  SM  (I) 
XB 1  I  )  3 X B I  I > 

YB I  I )=YB<  I )*SUNBET 
ZB  I  I ) =ZB 1  I )*SUN8ET 
SMI  I  ) =SM I  I ) 

102  CONTINUE 

IF  I  ISM  .GT.  0  >  GO  TO  104 
IM  3  12 

DO  103  1=1,  IM 
SK  I  I  )  3  SM  I  I  ) 

VK  (  I )  3  VS  I  I ) 

XK  (I  )  3  XB  I  I  ) 

ZK  I  I  >  3  ZB  II) 

103  CONTINUE 
GO  TO  150 

104  IF  (  12  -GT.  50  >  GO  TO  140 

DO  120  K  3  1 ,  IM 

SK  ( K )  3  SM  I  12 )  *  PCHORD  IK) 

120  CONTINUE 

IR  3 1 2  -  1 


CALL 

WAC  I 

12, 

SM. 

XB, 

IM.  SK. 

XK. 

OF,  TF 

CALL 

WAC  I 

12. 

SM. 

ZB  , 

IM,  SK, 

ZK, 

DF  .  TF 

CALL 

WAC  I 

12. 

SM, 

VB  . 

IM,  SK, 

VK, 

DF,  TF 

IFIKPRINT. 

.  EQ • 0 ) 

GO 

TO 

150 

WRITE  I  6,  51  ) 

DO  1 30  K  3  1 .  IM 
ZK(K)3ZK(K)/SUNBET 

WRITE  I  S,  52  >  K,  PCHORD  IK),  XK  <K>,  ZK  <K>,  SK  IK).  VK  (K) 
ZK(K)3ZK(K)*SUNBET 
130  CONTINUE 
GO  TO  150 
140  IM  3  12 

150  IF  I  IBETA  .EQ.  0  )  GO  TO  180 
BETAK  (1)  3  1.0 
I M 1  3  IM  -  1 

I  M2  3  IM  -  2 
DO  160  I  3  2,  I M 1 

D I  3  -I  SKI I-l ) -SKI  I ) >  /  <  ISKI  1*1 l-SKI  1-1 >  )  *  ( SK I  I >-SK( I  -  1 ) ) ) 

FI  3  ISKI  I  l-SKI  1-1 )  )  /  (  ISKI 1*1 )-SK< 1-1 ) )  *  I SK I  I  +  1 ) -SK I  I ) )  ! 

El  3  - 1  DI  ♦  F  I  ) 

DVDS  -  01  *  VK  I  I  - 1 )  »  E I  3  VK  (I)  *  F I  *  VK  fl  +  l) 

BETAK  I  I )  3  I  SK  I  I )  /  VK  I  I )  )  *  DVDS 
160  CONTINUE 

PI  -  SK  I IM)  -  SK  I IM1 ) 

P  2  3  SK  I  IM )  -  SK  I  I  M2 ) 

P  3  3  SK  I  IM1  )  -  SK  I IM2 ) 

DI  3P1/(P3*PZ) 

FI  3  P  2  /  I  P  3  *  P  1  ) 

El  3  I  2.0*SK(IM)  -  SKI IM1 )  -  SKIIMZ)  )  /  (  P 1  •  P2  ) 

DVDS  3  DI  *  VKIIM2)  -  FI  *  VKIIMl)  *  El  *  VKIIM) 

BETAK  (IM)  3  I  SK  (IM)  /  VK  (IM)  )  *  DVDS 
WRITE  I  6,  54  ) 

DO  165  I  3  I,  IM 
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o  o  o  o 


DO  210  I  1  1 ,  KSORCE 
ZC< I >=ZCt I ) /SUNSET 

WRITE  (  6,  41  )  LK,  I,  XC  <I),  ZC  t I > .  V  (I),  S  (I). 

I  D(I>.  ODDS  (  I  ) 

ZC( I )=ZC( I )»SUNBET 
2 10  CONTINUE 

IF  (  NOCAL  .EQ.  0  )  GO  TO  248 
240  DO  242  K  =  1  ,  KSORCE 
DODS  ( K )  =0.0 
D  (K)  *  0.0 
242  CONTINUE 

THIS  PARTICULAR  ASSIGNMENT  OF  KK=KK+KSORCE  IS  TO  INCORPORATE 
THE  ELEMENT  COUNT  WHILE  THE  NON-LIFTING  PART  IS  NOT  TAKEN 
FOR  VISCOUS  CALCULATIONS .THIS  IS  DISCARDED  HERE. 

KK  =  KK  ♦  KSORCE 
248  CONTINUE 

C  PUT  THE  COMPUTED  RHS  TERMS  IN  THE  ORDER  OF  THE  CONTROL  POINTS. 

C 

KR  =0 

DO  250  J  =  Ji ,  J2 
KR  =  KR  +  1 

RHS  (J)=  OODS  <  KR  ) 

250  CONTINUE 

KTOTAL  =  KTOTAL  +  KSORCE 
C 

C  INCREMENT  TO  THE  NEXT  STRIP 

C 

300  CONTINUE 

LI  =  L2  *  1 

GO  TO  1500 
C 

C  THE  FOLLOWING  PART  IS  FOR  THE  LIFTING  STRIPS. 

C 

500  I  ERR  =  0 

IM  =51 

KSTRIP  =  0 

KN  •  0 

NL  =0 

K  1  =  LI 

LP  =0 

L  S  *0 
L  S  =  L  S  +  1 

NS  »  NSTRIP  I  LS  ) 

I  EXTRA  *  I XF  LAG  (  LS  ) 

KWAKE  =  NWAKE  (  LS  ) 

IPOINT  =  NSORCE  (  LS  ) 

NP  =  KWAKE  +  IPOINT 

IF  (  I  EXTRA  .Ed.  0  >  GO  TO  530 

IF  (  I  EXTRA  .EQ.  1  .OR.  I  EXTRA  .EQ.  3  )  GO  TO  510 
K2  =  L2  ♦  2 

GO  TO  550 

5 10  K2  *  L2  ♦  1 

GO  TO  550 
530  K2  «  L  2 

550  DO  1 300  LK  =  LI,  L2 
KN  =  KN  +  1 

LP  *  LP  ♦  1 

IPOINT  *  NSORCE  <  LS  ) 

LIG1  »  I G 1  !  LP.  LS  > 

LIGN  =  IGN  (  LP.  LS  ) 

NIGNOR  =  LIGN  -  I.IG1 

IF  (  LIGN  .EQ.  0  .AND.  LIG1  .EQ.  0  >  GO  TO  600 
NIGNOR  =  NIGNOR  ♦  1 

600  KSORCE  *  IPOINT  -  NIGNOR 

Jl  =  KK  +  1 

J 2  »  KK  +  KSORCE 

VB  ( 1 )  =0. 

DELS ( 1 >=  0. 
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SM  (  1) 

- 

DELS  (1) 

DO  650 

K 

■  1.  KSORCE 

KK 

a 

KK  +  1 

A  1  1 

» 

UNITVX 

{ 

KK  ) 

A I  2 

* 

UNITVY 

{ 

KK  ) 

A  1  3 

a 

UNITVZ 

( 

KK  ) 

VX 

= 

VELX 

< 

KK  ) 

VY 

= 

VELY 

( 

KK  ) 

VZ 

» 

VELZ 

( 

KK  ) 

V  (K) 

* 

VEL 

( 

KK  ) 

XC  (K) 

* 

XCENT 

( 

KK  ) 

YC  (  K  ) 

* 

YCENT 

( 

KK  ) 

ZC  (K) 

* 

ZCENT 

( 

KK  ) 

T  (  K  ) 

» 

All  * 

VX 

+  A 1 2 

CONTINUE 

KM  1 

a 

KSORCE 

- 

1 

DO  700 

K 

=  1.  KM  1 

TPROD 

a 

T  (  K 

) 

*  T  ( 

VY  *  A1 3 


IF  (  TPROD  . GE.  ZERO  >  GO  TO  700 
KSTART  =  K 


vz 


GO  TO  720 
700  CONTINUE 


C 

C  NO  SEPARATION  OF  FLOWS,  WRITE  A  MESSAGE  AND  QUIT. 

C 

WRITE  (  6.  80  ) 

GO  TO  4000 


FOLLOWING  STATEMENTS  ARE  FOR  GOOD  RUNS. 
7Z0  CONTINUE 


DENOM 

-  T  ( 

KSTART  + 1  ) 

- 

T  ( 

KSTART  ) 

XI 

-  (  T 

<  KSTART* 1 

) 

*  XC 

(  KSTART  ) 

- 

1 

T 

(  KSTART  ) 

* 

XC  ( 

KSTART* 1  ) 

>  / 

DENOM 

YI 

*  (  T 

(  KSTART* 1 

) 

*  YC 

(  KSTART  ) 

- 

1 

T 

(  KSTART  ) 

* 

YC  ! 

KSTART*  1  > 

>  / 

DENOM 

ZI 

•  (  T 

<  KSTART* I 

) 

*  ZC 

(  KSTART  ) 

- 

1 

T 

(  KSTART  ) 

« 

ZC  < 

KSTART+1  ) 

>  / 

DENOM 

XI.YI.ZI  ARE  COORDINATES  OF  THE  INITIAL  POINT  ON  BOTH  LIFTING 
BOUNDARY  LAYER  STRIPS. 


VB  { 1 >  =0. 

DO  1200  IBL  -  1.  2 
READ  (  5.  97  )  NXT 
CALL  HEADER 
WRITE  (  6,  10  ) 

WRITE  (  6,  20  )  IS,  JTYPE ,  JSTRIP 
X!«XI 

Y I «Y I /SUNBET 
ZI-ZI/SUNBET 

WRITE  (  6.  70  )  LK.  KSTART ,  XI,  YI ,  ZI 

YI*YI*SUNBET 

ZI-ZI*SUNBET 

REWIND  13 

IB  »  IBL 

XB  ( 1 )  •  XI 

YB  ( 1 )  -  YI 

ZB  ( I >  -  ZI 

IF  (  18  .EQ.  2  >  GO  TO  800 
12  •  KSTART  +  1 

KM  *  1 

DO  750  K  -  2,  12 

KL  -  12  -  KM 

XB  ( K)  .-  XC  (  KL  > 

YB  (K)  -  YC  <  KL  ) 

ZB  <K>  ■  ZC  (  KL  > 

VB  (K>  »  V  (  KL  ) 

KM  -  KM  ♦  1 
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750  CONTINUE 
GO  TO  910 
800  CONTINUE 

12  -  KSORCE  -  KSTART  ♦  1 

KM  ■  KSTART 

OO  900  K  -  2.  12 

KL  -  KM  ♦  1 

XB  (K)  »  XC  (  KL  ) 

YB  (<)  ■  YC  (  KL  ) 

ZB  (K>  *  2C  (  KL  > 

VB  ( K  )  *  V  (  KL  ) 

KM  *  KM  ♦  1 

900  CONTINUE 

C 

C  CALL  CALC3L  TO  CALCULATE  ARC  LENGTHS  AND  CALL  BL .  PROGRAM. 

C 

910  DO  915  K  *  2.  12 

DELX  *  <  XB  (  K  )  -  XB  (  K-l  >  >  *•  2 

DELV  =  (  YB  (  K  >  -  YB  (  K-l  )  )  **  2 

DeLZ  =  (  Z8  <  K  1  -  ZB  <  K-l  )  >  **  2 

DELS(K>=  SORT  (  DELX  *  DELY  *  DELZ  > 

SM  <K)  =  SM  (  K-l  )  ♦  DELS  <  K  > 

915  CONTINUE 

WRITE  (  6.  30  ) 

DO  920  1=1.  12 
YB 1  I >»Y8<  I  l/SUNBET 
ZB ( I >=ZB< I l/SUNBET 

WR  I  TE  <  6 . 9  9  1  BLIND(IBL).  I.  XBII),  YB  <  I  )  ,  ZB(I).  VB  (  I  )  ,  SM(  I  > 
YB( I >*YB< I )*SUNBET 
ZB( I >=Z8( I >*SUNSET 
920  CONTINUE 

IF  (  ISM  .GT.  0  )  GO  TO  928 
IM  «  12 

DO  924  I  =  1  .  IM 
SK  (I)  =  SM  (  I  ) 

VK  (I)  =  VB  (  I  ) 

XK  (I  1  =  XB  (II 
ZK  ( I )  =  ZB  (  I  ) 

924  CONTINUE 
GO  TO  950 

928  IF  (  12  .GT.  50  )  GO  TO  940 
DO  930  K  *  1 ,  IM 
SK  (K)  =  SM  (12)  *  PCHORD  (K) 

930  CONTINUE 

IR  =12-1 


CALL 

WAC  ( 

12. 

SM. 

XB, 

IM, 

SK, 

XK, 

OF, 

TF 

CALL 

WAC  ( 

12  . 

SM. 

ZB  , 

IM, 

SK, 

ZK, 

OF  , 

TF 

CALL 

WAC  < 

12, 

SM, 

VB  . 

IM, 

SK, 

VK, 

DF  , 

TF 

I F ( KPR I  NT . EO . 0 )  GO  TO  950 
WRITE  (  6,  51  ) 

DO  932  K  -  1 ,  IM 
ZK(K)=ZK(K)/SUN8ET 

WRITE  (  6.  52  )  K,  PCHORD  (K),  XK  (K>,  ZK  (K),  SK  (K),  VK  (K) 
ZK(K)=ZK(K)*SUNBET 
932  CONTINUE 
GO  TO  950 

940  IM  =  12 

950  IF  (  IBETA  .EO.  0  >  GO  TO  970 
BETAK  ( 1  1  ■  1.0 

I M 1  •  IM  -  1 

I  M2  =  IM  -  2 

DO  955  I  =  2,  I M 1 

01  »  - ( SK( I  +  1 > -SK < I  ) )  /  ( (SKI  1*1 >-SK( 1-1  >  )  •  (SK(I)-SK(I-l))) 

FI  -  (SK( I )-SK( 1-1 1  1  /  ( (SK( 1  + I ) -SK ( 1-1)  )  *  < SK ( I ♦ I  ) -SK( I)  )  ) 

El  «  -<  DI  ♦  FI  > 

DVOS  »  DI  *  VK  (  I  - 1 )  ♦  El  «  VK  (I)  +  F I  »  VK  (1*1) 


i 
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WRITE  (  6.  55  )  I.  SK  (I),  VK  <I>,  BETAK  (I) 

165  CONTINUE 
180  CONTINUE 

WRITE  <  13  )  IM,  IBETA 

IF  <  IBETA  .GT.  0  )  WRITE  (  13  )  (  BETAK  (K).  K  «  1.  IM  > 


WRITE  ( 

13 

(  XK 

(K)  , 

K  -  1 

.  IM 

> 

WRITE  ( 

13 

(  ZK 

<  K  )  , 

K  -  1 

.  IM 

> 

WRITE  < 

13 

(  VK 

(K)  . 

K  -  1 

.  IM 

) 

WRITE  ( 

13 

(  SK 

(K)  , 

K  -  1 

.  IM 

) 

WRITE  ( 

1 3  : 

NXT 

REWIND  13 
C 

IF(KPRINT.GT.0>  WRITE  <  6 .  96  > 

CALL  80UNDL 
REWIND  2 
READ  (  2  )  IK 

READ  <  2  >  (  DEIK  <K>.  K  -  1.  IK  > 

REWIND  2 

IF  (  IK  .EQ.  IM  )  GO  TO  182 
WRITE  <  6,  94  )  IM,  IK 
STOP 

182  CONTINUE 

IF  (  ISM  .EQ.  0  )  GO  TO  184 
C 

CALL  WAC  (  IK,  SK,  DEL K ,  12,  SM.  DEL,  DF ,  TF  ) 

GO  TO  186 

184  DO  185  K  *  1 ,  IM 
DEL  (K)  •  DELK  (K) 

185  CONTINUE 

IF1KPRINT.EO.0)  GO  TO  189 

186  WRITE  (  6,  53  ) 

DO  188  K  =  1  .  12 

WRITE  (  6,  55  )  K.  SM  !K>,  VB  <K>,  DEL  (K) 

188  CONTINUE 
189  CONTINUE 

C 

IR-KSORCE 
CALL  CALCRS 

DO  200  K  -  2,  12 
S  (K-l )-  SM  (  K  ) 

D  (K-l  >-  DEL  <  K  ) 

V  (K-l >-  VB(  K  ) 

DQDS  (K-l )  -  DOS  (  K  > 

200  CONTINUE 

GO  TO  ( 1 1 120, 1 1 130) , IBB 

11120  KM-KBSTRT 
DO  11121  K-2, 12 
S ( KM ) -SM ( K  ) 

0 ( KM ) -DEL ( K ) 

V ( KM ) -VB ( K  > 

DQDS(KM)-OQS(K) 

KM-KM-1 

11121  CONTINUE 

GO  TO  11115 

11130  KM-KBSTRT-1 

DO  11122  K-2,  12 
S ( KM  >  *SM( K ) 

D  (  KM ) -DEL ( K  > 

V ( KM ) -VB ( K ) 

DQDS ( KM ) -DOS ( K ) 

KM-KM-1 

11122  CONTINUE 

11115  CONTINUE 

CALL  HEADER 

IF(KPRINT.GT,0)  WRITE  (  6 ,  50  ) 

IF(KPRINT.GT.0>  WRITE  (  6 ,  60  >  (  DQDS  Cl),  I  -  I,  KSORCE  ) 
WRITE  (  6,  32  ) 


129 


I 


BETAK  I  I  >  -  I  SK  I  I >  /  VK  I  1 >  >  *  DVDS 
955  CONTINUE 

PI  ■  SK  { IM>  -  SK  ! IM1 ) 

P 2  -  SK  I IM)  -  SK  < 1  M2 > 

P 3  =  SK  (  IM1 )  -  SK  <  IM2) 

D I  -  PI  /  (  P3  *  P2  ) 

FI  *  P2  /  I  P3  *  P 1  > 

El  *  (  2.0-SKIIM)  -  SKI  I  HI >  -  SKIIM2)  >  /  (  P1*P2  ) 

DVDS  -  01  *  VK I  I  M2 )  -  FI  *  VK(IMl)  ♦  El  *  VK(IM) 

BETAK  I  I M  >  «  I  SK  <  I M  >  /  VK  <  I M  >  >  *  DVDS 
WRITE  (  6.  54  ) 

DO  965  I  -  1  .  I M 

WRITE  I  6,  55  )  I.  SK  (I),  VK  (I).  BETAK  (I) 

965  CONTINUE 
970  CONTINUE 
C 

WRITE  I  13  )  IM,  IBETA 

IF  (  IBETA  .GT.  0  )  WRITE  {  13  )  (  BETAK  <K>,  K  *  1,  IM  > 


WRITE  I 

13  ! 

I  XK 

IK), 

K  -  1 

.  IM 

) 

WRITE  ( 

13  ; 

I  ZK 

IK), 

K  »  1 

,  IM 

> 

WRITE  I 

13  ; 

I  VK 

IK), 

K  =  1 

,  IM 

> 

WRITE  < 

1 3  ; 

I  SK 

IK), 

K  -  1 

,  IM 

> 

WRITE  < 

1 3  i 

NX  T 

REWIND  13 
C 

IFIKPRINT.GT.0)  WRITE  I  6,  96  ) 

CALL  BOUNDL 
REWIND  2 
READ  (  2  )  IK 

READ  (2)1  DELK  (K),  K  -  1.  IK  ) 

REWIND  2 

IF  I  IK  .EQ.  IM  )  GO  TO  1000 
WRITE  I  6,  94  )  IM,  IK 
STOP 

1000  CONTINUE 

IF  I  ISM  .EO.  0  )  GO  TO  1010 
C 

CALL  WAC  (  IK,  SK,  DELK,  12,  SM,  DEL,  DF ,  TF  > 
GO  TO  1020 


C 

1010  DO  1015  K  »  1 ,  IM 
DEL  IK)  *  DELK  (K) 

1015  CONTINUE 
1020  CONTINUE 

IFIKPRINT.EQ.0)  GO  TO  1026 
WRITE  (  6,  53  ) 

DO  1025  K  »  1.  12 

WRITE  (  6,  55  )  K,  SM  IK),  VB  IK),  DEL  IK) 

1025  CONTINUE 

1026  CONTINUE 


C 

CALL  CALCRS 
C 

IR  *  KSORCE 

C 

C  PUT  THE  S"S,  DEL "S  BACK  TO  ORDER. 
C 

GO  TO  I  1050,  1 100  ) ,  IB 
1050  KM  ■  KSTART 

DO  1060  K  •  2,  12 


S  (KM)  -  SM  I  K  ) 


1 
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0  ( KM  >  »  DEL  (  K  ) 

V  (KM)  -  VS  (  K  ) 

DQDS  (KM)  -  DOS  (  K  ) 
KM  -  KM  -  1 

1060  CONTINUE 

GO  TO  1200 

1100  KM  >  KSTART  ♦  1 

DO  1 150  K  -  2,  12 
S  (KM)  •  SM  (  K  > 

D  (KM)  *  DEL  (  K  ) 

V  (KM)  *  VB  (  K  > 

DQDS  (KM)  -  DOS  (  K  ) 
KM  •  KM  +  1 

1150  CONTINUE 
1200  CONTINUE 


CALL  HEADER 

I F ( KPR I  NT . GT .0 )  WRITE  (  6 .  50  ) 

I F ( KPR I  NT . GT . 0 )  WRITE  !  6 ,  60  )  <  DQDS  <I),  I  »  1,  KSORCE  ) 
WRITE  (  6.  32  ) 

DO  1220  1-1,  KSORCE 
ZC( I )-ZC( I ) /SUNSET 

WRITE  (  6.  41  )  LK,  I,  XC  (I),  ZC  (I),  V  (I),  S  (I), 

1  D  (  I ) ,  DQDS  ( I  ) 

ZC ( I >  *ZC ( I ) “SUNSET 
1220  CONTINUE 

PLACE  THE  COMPUTED  RHS*S  IN  THE  ORDER  OF  THE  CONTROL  POINTS. 

KR  *  0 

DO  1250  0  -  01 ,  02 
KR  »  KR  +  1 

RHS  (0)-  DQDS  (  KR  ) 

1250  CONTINUE 

KTOTAL  -  KTOTAL  ♦  KSORCE 
1300  CONTINUE 

LI  *  L2  ♦  1 

1500  CONTINUE 

IF  (  KTOTAL  .EQ.  KONTRL  )  GO  TO  2000 
WRITE  (  6.  92  ) 

GO  TO  4000 

PUT  THE  NEW  RHS  COLUMN  ON  UNIT  IRSIDE. 

2000  WRITE  (  IRSIDE  )  IONE 

WRITE  (  IRSIDE  >  (  RHS  (  I  ),  I  »  1,  KONTRL  ) 

WRITE  (  $,  50  ) 

WRITE  <  6,  60  )  <  RHS  (  I  ).  I  *  1.  KONTRL  ) 

REWIND  IRSIDE 
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CSC*CECK  PRINT 

OVERLAY (LINK, 6,2) 

L 1ST, NONE 
PROGRAM  PR  T NT 
SUBROUTINE  PRINT 

THIS  SUBROUTINE  PRINTS  OUT  THE  FINAL  OUTPUT  WHICH  HAS  THE  MOST 
IMPORTANT  DATA. 


LOGICAL 

INTIAL 

LAST  , 

SKIP 

REAL 

MOMENX 

MOMENY, 

MOMENZ 

REAL 

MX  I 

( 

500 

)  , 

MYU  ( 

500 

) 

REAL 

MZK 

! 

500 

) 

DIMENSION 

COMSIG 

( 

1  100 

>  . 

AREA  < 

1 100 

1 

DIMENSION 

VX 

I 

1100 

)  , 

VY  ( 

1 100 

) 

DIMENSION 

VZ 

( 

1  100 

) . 

VSUM  ( 

1  100 

) 

DIMENSION 

VEL 

< 

1  100 

)  , 

OCX  < 

2  100 

) 

DIMENSION 

DCY 

I 

1100 

)  . 

DCZ  ( 

1  100 

) 

DIMENSION 

VN 

{ 

1  100 

)  . 

PCOEF  ( 

1100 

> 

DIMENSION 

NSORCE 

I 

10 

)  , 

NWAKE  I 

10 

) 

DIMENSION 

IG1  (  50 

.  10 

> , 

IGN  (  50, 

10 

) 

DIMENSION 

NLINE 

( 

50 

)  , 

NSTRIP  ( 

10 

) 

DIMENSION 

NLT 

( 

500 

>  , 

NTYPE  ( 

50 

) 

DIMENSION 

NL INE1 

< 

10 

)  , 

NLINEN  ( 

10 

) 

DIMENSION 

KSPL  IT 

I 

10 

)  . 

KOUNT  ( 

50 

) 

DIMENSION 

FSTRPX 

I 

500 

)  , 

FSTRPY  ( 

500 

) 

DIMENSION 

FSTRPZ 

( 

500 

)  , 

FSECX  ( 

50 

> 

DIMENSION 

FSECY 

( 

50 

>  , 

FSECZ  ( 

50 

) 

DIMENSION 

SECMOX 

( 

50 

> , 

SECMOY  ( 

50 

) 

DIMENSION 

XC 

( 

1100 

) 

,  VC 

t 

1  100 

DIMENSION 

zc 

I 

1100 

) 

XN 

( 

1100 

DIMENSION 

YN 

< 

1 100 

) 

,  ZN 

( 

1100 

DIMENSION 

SECMOZ 

< 

50 

> 

DIMENSION 

ALPHAX 

l 

10  ) 

f 

ALP HAY  ( 

10  > 

DIMENSION 

ALPHAZ 

( 

10 

)  , 

TITLE  ( 

18 

) 

COMMON 

/ 

CONFLG 

/ 

TITLE, 

CASE,  LIFSEC.  LIST,  NOFF. 

I OUT ,  KONTRL 

KUTTA, 

MPR 

COMMON 

/ 

INFORM 

/ 

NSORCE. 

,  NWAKE.  NSTRIP, 

I G 1 ,  IGN 

COMMON 

/ 

BOINFO 

/ 

KOUNT. 

NLINE,  NTYPE,  NLT,  ISECT 

COMMON 

/ 

IMPORT 

/ 

NON,  NF  LOW 

COMMON 

/ 

VBLOCK 

/ 

VX,  VY, 

VZ,  VEL 

COMMON 

/ 

FINAL 

/ 

COMSIG, 

VSUM,  DCX,  DCY, 

DCZ,  VN. 

PCOEF,  B ( 50> 

COMMON 

/ 

EXTRA 

/ 

IXFLAG 

(  10  ),  SKIP 

COMMON 

/ 

ORIGIN 

/ 

ORIGNX. 

ORIGNY,  ORIGNZ 

COMMON 

/ 

ANGLE 

/ 

MANGLE . 

ALPHAX,  ALPHAY, 

ALPHAZ 

.SUNBET 

COMMON 

/ 

ATTACK 

/ 

ANGLEX . 

ANGLEY,  ANGLEZ 

COMMON 

/ 

OUTPUT 

/ 

F  PR  I  NT , 

IANGLE 

COMMON 

/ 

CTABLE 

/ 

XC,  YC, 

ZC  .AMACH 

COMMON 

/ 

NORMAL 

/ 

XN,  YN 

ZN 

COMMON 

/ 

ATABLE 

/ 

AREA 

COMMON 

/ 

SCALE 

/ 

FC,  BL 
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100  F OR MAT ( 1  HO , 50X , 32HON  -  BODY  POINTS  FINAL  OUTPUT  / 

1  1H  ,50X  .  32  C  1H-)  ) 

110  F  OR MAT ( 1 H  .11X.4HN  M , SX , 2HX0 , 8X , 2HY0 , 8X . 2HZ0 , 8X , 2HVX , 

1  8X , 2HVY , 8X , 2HVZ , 8X , 2HC  P , 8X , 2H VN , 5X , 5  HS I GMA ) 

120  F  OR MAT ( 1 H  , 9X , 2 1 3 . 9F 10 . 6 > 

130  F  OR MAT ( 1H  , 1 2X , 1 3 . 9F 10 . 6 > 

140  FORMAT ( 1 M  , 8X , 3 ( 20X . 1 2H *«*»«*••****> > 

150  FORMATUH  ,  29X  ,  20HSTR  I P  (  FORCE/Q  )  ,  3 (  2X  ,  E  1 3  .  S  >  > 

155  FORMATUH  ,  29X  ,  20HSECTION  <  FORCE/Q  )  .  3  (  2X  .  E  1  3 . 6  )  > 

158  FORMATUH  .29X.20HBODY  (  FORCE/Q  >  ,  3  (  2X  .  E 1  3 . 6  >) 

160  F  ORMAT ( 1 H  , 29X , 20HSTR I P  (  MOMENT /Q  ).3(2X.£13.6>) 

162  FORMATUH  .  29X  ,  20HSECT I  ON  (  MOMENT/Q  >  ,  3  (  2X  .  E  1 3 . 6  >  > 

164  FORMATUH  .29X.20HBODY  (  MOMENT/Q  )  ,  3  (  2X  ,  E  1  3 . 6  )  ) 

C 

200  FORMAT  (  1H0,  52X .  28HKUTTA  POINTS  FINAL  OUTPUT  / 

1  1H  .  52X ,  28  (  1H-  )  ) 

210  FORMAT  (  1H0,  32X .  5HKUTTA,  7X ,  2HX0,  11X.  2HVX ,  11X.  2HVT ,  11X, 

1  3HDCX .  11X,  2HNX  /  1H  ,  32X,  6HPOINTS,  6X.  2HY0,  11X, 

2  2HVY,  11X,  2HVN.  11X,  3HDCY,  UX.  2HNY  /  1H  .  44X,  2H20. 

3  11X.  2HVZ,  11X,  2HCP .  1 1 X .  3HDCZ,  11X,  2HNZ  /  1H  .  32X, 

4  6  <  1 H-  )  ,  3X,  9UH->.  4(  4X  .  9  (  1 H  -  >  >  ) 

220  FORMAT  (  1H0,  33X.  14,  5F13.6  /  (  38X.  5F13.6  1  ) 

230  FORMAT  (  1H0,  51X.  9HSTRIP  NO..  10X,  11HB  (  STRIP  >  / 

1  1H  .  SIX.  9!  1H-)  .  9X ,  1 3 (  1 H-  )  /  ) 

240  FORMAT  <  1H  .  53X .  14.  12X,  E13.6  ) 

300  FORMAT  (  1H0.  49X.  33HOFF  -  BODY  POINTS  FINAL  OUTPUT  / 

1  1H  ,  50X.  33 (  1H-  )  ) 

310  FORMAT  (  1H0,  38X .  6HPOINTS,  7X .  2HX0,  11X,  2HVX ,  12X.  2HVT. 

1  10X.  3HDCX.9X. *  UOFOUT '  /  1H  .  51X.  2HY0.  11X.  2HVY.  UX ,  4HVTSQ , 

2  9X ,  3HDCY.9X, 'VOFOUT'  /  1H  ,  51X,  2HZ0.  1 1 X .  2HVZ,  1 2X .  2HCP , 

3  10X,  3HDCZ.9X,  'VOFOUT*  /  1H  ,  38X.  6I1H-),  5<  4X ,  9UH-) 

4  )  ) 

320  FORMAT  (  1H0.  39X .  13.  2X,  SF 1 3 . 6  /  t  45X.  5F13.6  )  > 

330  FORMAT  (  1 H0 ,  40X .  I4.28H  ON  -  BOOY  POINTS  MISSED.  . 

1  22HEXECUTION  TERMINATED.  ) 

340  FORMAT  (  1 H0 ,  41X,  14,  24H  KUTTA  POINTS  MISSED.  . 

1  22HEXECUTION  TERMINATED.  ) 

350  FORMAT  (  1H0.  40X .  14,  27H  OFF-BODY  POINTS  MISSED.  . 

1  22HEXECUTION  TERMINATED.  > 

C 

C  -  * - - 

c 

OUTPUT  INCLUDES  3  PARTS-- 
C  (1)  THE  ON-BODY  POINTS  FINAL  OUTPU1', 

C  <2'  THE  KUTTA  POINTS  FINAL  OUTPUT.  AND 

C  (3)  THE  OFF-BODY  POINTS  FINAL  OUTPUT 

r 


Z  KT25W  -  PRESSURE  DATA  OUTPUT  FLAG  <0  TO  SUPPRESS) 

C 

KT25W- 1 

SKIP  -  .FALSE. 

C 

CALL  HEADER 
VRITE16. 100) 

WR ITE ( 6 . 1 10) 

NPRT  -  3 
C 

L  1  •  1 

LS  -  0 
K  K  ■  0 

KM  1  •  NFLOV  -  MANGLE 
KN  •  0 
FORCEX  ■  0. 

FORCEY  -  0. 

FORCEZ  -  0. 

MOMENX  •  0. 

MOMENY  *  0. 

MOMENZ  ■  0. 
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FC2 

• 

FC  * 

FC 

DO  360 

NY 

1  -  I 

,  500 

MXr  ( 

NY 

)  - 

0. 

MY  J  ( 

NY 

)  - 

0. 

MZX  ! 

NY 

)  - 

0. 

FSTRPX 

( 

NY  ) 

•  0. 

FSTRPY 

( 

NY  ) 

*  0. 

FSTRPZ 

( 

NY  ) 

-  0. 

DO  362 

NR 

«  1 

.  50 

SECMOX 

( 

NR  ) 

-  0. 

SECMOY 

( 

NR  ) 

-  0. 

SECMOZ 

! 

NR  ) 

*  0. 

F  StCX 

!  NR  ) 

-  0. 

FSECY 

( 

NR  ) 

*  0. 

FSECZ 

{ 

NR  ) 

-  0. 

NPRT  - 

-  PRINT 

-LINE 

»**.*L1  -  INITIAL  STRIP  COUNT 

-  LIFTING  SECTION  INCREMENT  FLAG 

*****KK  -  CONTROL  ELEMENT  INCREMENT  COUNT 

PPPP-0. 

DO  1000  IS  -  1,  ISECT 
J 1SUND-NTYPE! IS ) 

J2SUN0-NL I NE ( IS) 

J  3SUND -XOUNT !  IS  ) 

WRITE!6.9923) 

9923  FORMAT  <  2X , 'THE  CONFIGURATION  I S  AF DL . 24 , MACH-0 . 925 • ) 

CS  PUNCH  9923 

VRITE<7,9923) 

WR ITE ! 6 , 5000)  IS . J1SUN0 , J3SUND, J2SUND 
C  PUNCH  5002, IS.J1SUN0.03SUND, J2SUND 

WRITE! 7,5002)  IS.J1SUN0.J3SUND.D2SUND 
FORMAT! 1H0, 1 IX. 'BODY  SECTION  NO- 
1  10X,;TOTAL  NO  OF  POINTS-’.  I4.10X 


5000 

5002 


1 2 , 1 0X  .  ’  TYP  E  -  ’ 
NO.  OF  STRIPS 


I  1  , 


_  ’Z"  '  wr  rwnng-  «  4«*,i»/\,  nw,  Uf  i  I  K  1  P  S  ■’  13  ) 

FORMAT! 2X, ’BODY  SECTION  NO- ’  , 13 . 3X . 'TYPE- ’  . I  3 . 3X , ’ TOTAL  NUMBER  0 
F  POINTS- ’, 15 , 3X ,’ NUM8ER  OF  STRIPS-’. 13) 


364 

366 

368 

369 

370 

380 

385 


1 

JTYPE  -  NTYPE  !  IS  ) 

IF  !  JTYPE  .EQ.  0  )  GO  TO  380 
LP  -  0 
LS  -  LS  *  1 

IPOINT  •  NSORCE  (  LS  ) 

NXTRA  -  IXFLAG  (  LS  ) 

IF  (  NXTRA  .EQ.  0  )  GO  TO  369 
IF  (  NXTRA  -  2  )  364.  366,  368 
INTIAL  -  .TRUE. 

LAST  -  .FALSE. 

GO  TO  370 
INTIAL  -  .TRUE. 

LAST  -  .TRUE. 

GO  TO  370 
INTIAL  -  .FALSE. 

LAST  -  .TRUE. 

GO  TO  370 
INTIAL  -  .FALSE. 

LAST  -  .FALSE. 

OSTRIP  •  NSTRIP  (  LS  ) 

GO  TO  385 

OSTRIP  •  NLINE  <  IS  ) 

12  -  LI  ♦  OSTRIP  -  1 
DO  900  LX  -  LI .  L 2 
KN  -  XN  ♦  1 
IF  (  JTYPE  .EQ 
IF  !  (  LX  .EQ.  LI 
IF  (  (  LX  .EQ.  L  2 
LP  -  LP  +  1 
LIGI  -  IG1  (  LP.  LS  ) 

LIGN  -  IGN  (  LP,  LS  > 


0  )  GO  TO  420 


•AND.  INTIAL  )  GO  TO  900 
■AND.  LAST  )  GO  TO  900 


NIGNOR 


LIGN 


LIGI 
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IF  (  LIG1  .EQ.  0  .AND.  LIGN  ,EQ.  0  I  GO  TO  400 
NIGNOR  ■=  NIGNOR  +  1 
400  KSORCE  =  IPOINT  -  NIGNOR 
GO  TO  500 

420  KSORCE  =  IASS  <  NLT  {  KN  )  ) 

500  DO  800  K  =  1  ,  KSORCE 
KK  -  KK  +  1 
X0  *  XC  (  KK  ) 

V0  -  YC  (  KK  ) 

Z0  *  ZC  (  KK  ) 

XNORM  -  XN  (  KK  ) 

YNORM  -  YN  (  KK  ) 

ZNORM  -  ZN  (  KK  ) 

C 

C  CALCULATE  THE  FORCE  COMPONENTS. 

C 

CPAREA  *  -PCOEF  < KK >  *  AREA  ( KK > 

FX  =  CPAREA  *  XNORM  *  FC2 

FY  =  CPAREA  «  YNORM  «  FC2 

FZ  =  CPAREA  *  ZNORM  *  FC2 

FSTRPX  (  LK  )  =  FSTRPX  (  LK  )  +  FX 

FSTRPY  (  LK  )  =  FSTRPY  (  LK  )  +  FY 

FSTRPZ  (  LK  )  =  FSTRPZ  <  LK  >  ♦  FZ 

C 

C  CALCULATE  THE  MOMENT  COMPONENTS. 

C 


RX 

X0  * 

FC  - 

ORIGNX 

RY 

= 

Y0  * 

FC  - 

ORIGNY 

RZ 

= 

Z0  * 

FC  - 

ORIGNZ 

XMO 

-  RY 

*  FZ 

-  RZ 

*  FY 

YMO 

-  RZ 

*  FX 

-  RX 

*  FZ 

ZMO 

-  RX 

*  FY 

-  RY 

*  FX 

MX  I 

(  LK 

1  »  1 

MX  I  ( 

LK  )  + 

XMO 

MYJ 

(  LK 

)  * 

MYJ  ( 

LK  )  + 

YMO 

MZK 

<  LK 

)  »  1 

MZK  ( 

LK  >  + 

ZMO 

XO 

a 

X0  * 

FC 

YO 

a 

Y0  * 

FC  /SUNBET 

ZO 

a 

Z0  * 

FC  /SUNBET 

AR 

a 

AREA 

(  KK 

)  *  FC 

*  FC 

6001  F  OR MAT ( 7F 10 . 5 , 12) 

IF  (  NPRT  .GE.  39  )  GO  TO  600 
NPRT  -  NPRT  ♦  1 
IF  <  K  .EQ.  1  >  GO  TO  700 

SUND 1  - VX <  KK ) 

SUND2*VY 1 KK ) 

SUN03=VZ!KK> 

CALL  CPSUNO(X0, Y0,Z0,SUND1 . SUN02 , SUND3 . SUNX , SUNY . SUNZ > 
WRITE  (  6.  130  )  K.  XO. YO.ZO. SUNX, SUNY, SUNZ, PCOEF(KK) , 

I  VN(KK) .COMSIG  (KK) 

C 

I F ( KT25W. EQ . 0 )  GO  TO  800 
KKW-0 
C 

GO  TO  800 
600  NPRT  -  0 

CS  CALL  HEADER 

CS  WRITE  (  6.  110  ) 

700  CONTINUE 

WRITE  (  6,  120  )  LK,  K.  XO , YO , ZO , SUNX , SUNY , SUNZ , PCOEF ( KK ) , 
1  VN(KK) .COMSIG  ( KK ) 

C 

I F ( KT25W . EQ . 0 )  GO  TO  800 
KKW-0 

IF(K.EQ.l)  KKW-I 
IFIK.EO.l  .AND.  LK.EQ.l)  KKV-4 
C 

800  CONTINUE 

IF  (  NPRT  .LT.  36  )  GO  TO  820 
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NPRT  -  0 
CALL  HEADER 
WRITE  (  6.  110  ) 

820  NPRT  -  NPRT  ♦  4 
WRITE  (  6.  140  ) 

WRITE  (  6.  150  )  PSTRPX  ILK),  FSTRPY  ILK),  FSTRPZ  ILK) 
WRITE  {  S.  160  >  MX I  (  LK).  MY  J  I  LK  >.  MZK  I  LK  ) 
WRITE  (  S.  140  ) 

FSECX  (IS)  -  FSECX  (IS)  ♦  FSTRPX  (LK) 

FSECY  (IS)  *  FSECY  (IS)  +  FSTRPY  (LK) 

FSECZ  (IS)  -  FSECZ  (IS)  ♦  FSTRPZ  (LK) 

SECMOX  (IS)  -  SECMOX  (IS)  +  MX  I  <  L  K ) 

SECMOY  (IS)  -  SECMOY  (IS)  ♦  MYJ  (LK) 

SECMOZ  (IS)  -  SECMOZ  (IS)  -  MZK  (LK) 

900  CONTINUE 

L  I  *  L2  *  1 

IF  <  NPRT  . LT .  37  )  GO  TO  920 
NPRT  =  0 
CALL  HEADER 
WRITE  (  6.  110  ) 

920  NPRT  »  NPRT  ♦  4 
WRITE  (  6.  140  ) 

WRITE  (  6.  1S5  )  FSECX  (IS),  FSECY  (IS).  FSECZ  (IS) 
WRITE  (  6.  162  )  SECMOX  (IS).  SECMOY  (IS),  SECMOZ  (IS) 
WRITE  (  6.  140  ) 

FORCEX  *  FORCEX  ♦  FSECX  (IS) 

FORCEY  -  FORCEY  *  FSECY  (IS) 

FORCEZ  -  FORCEZ  ♦  FSECZ  (IS) 

MOMENX  •  MOMENX  ♦  SECMOX  (IS) 

MOMENY  -  MOMENY  ♦  SECMOY  (IS) 

MOMENZ  -  MOMEN2  ♦  SECMOZ  (IS) 

1000  CONTINUE 

IF  (  NPRT  .LT.  37  )  GO  TO  1010 
NPRT  -  0 
CALL  HEADER 
1010  NPRT  -  NPRT  ♦  4 
WR  I TE  (  6 . 1  40 ) 

WRITE  (  6,  158  )  FORCEX,  FORCEY,  FORCEZ 
WRITE  <  6,  164  )  MOMENX,  MOMENY,  MOMENZ 
WRITE  <  6,  140  ) 

IF  (  KK  .NE.  KONTRL  I  GO  TO  3200 
IF  (  KUTTA  -EQ.  0  )  GO  TO  2100 

NOW  THE  KUTTA  POINTS  OUTPUT 


NPRT  *  45 

DO  2000  K-  1,  KUTTA 
KK  -  KK  ♦  1 
X0  -  XC  (  KK  ) 

Y0  ■  YC  (  KK  ) 

Z0  -  ZC  (  KK  ) 

XNORM  »  XN  (  KK  > 

YNORM  -  YN  (  KK  ) 

ZNORM  -  ZN  (  KK  ) 

IF  (  NPRT  .GE.  10  )  GO  TO  1500 
NPRT  «  NPRT  ♦  1 
GO  TO  1800 
1500  NPRT  -  0 

CALL  HEADER 
WRITE  (  6,  200  > 

WRITE  (  6.  210  > 


1 


X0«X0«FC 

Y0*Y0»FC/SUNBET 


Z0-Z0*FC/SUN8ET 

800  WRITE  (  6,  220  )  K,  X0,  VX !  KK  ),  VEL(  KK  ).  0CX<  KK  ),  XNORM 
,  Y0 ,  VY (  KK  ) ,  VN<  KK  ) ,  0CY(  KK  >.  VNORM,  Z0 

2  VZ(  KK  >.  PCOEFI  KK  ).  DCZ(  KK  ).  ZNORM 


X0-X0/FC 

Y0-Y0/FC*SUN8ET 
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Z0-Z0/FC*SUNBET 
2000  CONTINUE 

WRITE  (  6.  140  ) 

'NEXT,  THE  OFF-BODY  POINTS  OUTPUT 

2100  IF  (  NOFF  ,EQ.  0  )  GO  TO  3100 
CS  PUNCH  5003 

WRITE! 7, 5003) 

5003  FORMAT(5X, 'THE  FOLLOWING  PUNCHED  CARDS  ARE  FOR  THE  OFF  BODY  POIN 

I  TS  '  ) 

NPRT  -  45 

KCOUNT  «  KK  -  KONTRL 

IF  (  KCOUNT  ,NE.  KUTTA  )  GO  TO  3400 
OO  3000  K  *  l ,  NOFF 
KK  =■  KK  ♦  1 
X0  =  XC  (  KK  > 

V0  =  YC  (  KK  ) 

Z0  =  ZC  (  KK  ) 

IF  (  NPRT  .GE.  10  )  GO  TO  2500 
NPRT  =  NPRT  *  1 
GO  TO  2800 
2500  NPRT  »  0 

CALL  HEADER 
WRITE  <  6.  300  ) 

WRITE  (  6,  310  ) 

2300  SUN04  *VX { KK 1 

SUND5-VY! KK ) 

SUND6-VZ1 KK ) 

CALL  CPOFFIX0, Y0.Z0.SU NO  4, SUND5, SUNOS, UOFOUT.VOF  OUT , WOF  OUT ) 

X0*X0*FC 

Y0-Y0*FC/SUNBET 

Z0=Z0*FC/SUNBET 

WRITE  (  6,  320  )  K.  X0,  VX  f  KK  > .  V£ L<  KK  ),  OCX (  KK  I.UOFOUT. Y0, 

1  VY (  KK  ).  VSUM!  KK  ),  DCY(  KK  J.VOFOUT,  Z0.  VZ(  KK  ). 

2  PCOEF(  KK  ),  DCZ 1  KK  ) .WOFOUT 
X0-X0/FC 

Y0-Y0/FC*SUNBET 
Z0“Z0/FC*SUNBET 
3000  CONTINUE 

WRITE  (  6.  140  ) 

KCHECK  •  KK  -  KONTRL  -  KUTTA 
IF  (  KCHECK  ,NE.  NOFF  )  GO  TO  3600 
3100  CALL  HEADER 

WRITE  (  6,  230  ) 

DO  3110  K  *  1  .  KM  1 


WRITE 

<  6  , 

240  )  K,  B 

(  K) 

31  10 

CONTINUE 

WR  ITE 

t  6 . 

140  ) 

GO  TO 

4000 

3200 

MISS  - 

IABS 

(  KK  -  KONTRL  > 

WRITE 

(  6 . 

330  )  MISS 

GO  TO 

4000 

3400 

MISS  • 

IABS 

(  KCOUNT  - 

KUTTA  ) 

WR  ITE 

(  6  . 

340  1  MISS 

GO  TO 

4000 

3600 

MISS  * 

IABS 

(  KCHECK  - 

NOFF  ) 

WRITE 

(  6 . 

350  )  MISS 

4000 

RETURN 

END 
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INITIAL  DISTRIBUTION 


DT IC-DDAC 

AUL/LSE 

ASD/ENSZ 

AFATL/DLODL 

AFATL/CC 

HQ  USAF/SAMI 

OO-ALC/MHWRB 

AFIS/INT 

HQ  TAC/DRA 

HQ  USAFE/DOQ 

HQ  PACAF/DOQQ 

HQ  TAC/INAT 

ASD/XRX 

USA  TRADOC  SYS  ANAL  ACTY 
COMIPAC  (PT-2) 

HQ  PACAF/OA 

USA  BALLISTIC  RESCH  LAB 

AFATL/CCN 

AFATL/DLODA 

ASD/ENESS 

HQ  TAC/XPS  (STINFO) 

6510  TEST  GP/ENML 

USMC /RED STONE  SCI  INFO  CNTR 

NSC  (TECH  LIBRARY) 

SANDIA  NATIONAL  LABS 

NAVAL  WPNS  CNTR  (CODE  3243) 

AFATL/DLCA 

DEPT  OF  AEROSPACE  ENG 
INDIAN  INST  OF  SCIENCE 
UNIV  OF  TENNESSEE  SPACE 
INSTITUTE/DR  SUNDARAM 
UNIV  OF  TENNESSEE  SPACE 
INSTITUTE/TECH  LIB 
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